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ABSTRACT 


Dynamic Analysis of Arbitrary Profile Liquid Annular Seals 
and IVansient Analysis with Large Eccentric Motion. 


One of the mam objective* of this work is to develop a new dynamic anal- 
ysis for liquid annular seal, with arbitrary profile and analyne a general distorted 
interstage seal of the Space Shuttle Main Engine High Pressure Oxygen Turbopump 
(SSME-ATD-HPOTP). The dynamic analysis developed is based on a method orig- 
inally proposed by Nelson and Nguyen (1988a, 1988b). The original method used 
an approximation scheme based on Post Fourier Tronoform, (FFT) to compute the 
arcumferentiai gradients of the seroth order variables of the eccentric solution. This 
method, in some cases, has problems with convergence at higher eccentricities. A sim- 
pler schmne based on cubic splines is found to be compntationaDy more efficient and 
la. better convergence properties at higher eccentricities than the original method. 
The first order solution of the original analysis is modified by including a more exact 
solution that takes into account the variation of perturbed variables along the cir- 
enmference. A new set of equation, for dynamic analysis are derived based on this 
more general model. The original method was developed for Moody’s friction model. 

In the current work, a unified solution procedure that is valid for both Moody’s and 
Hirs’ models is presented. 

Dynamic analysis based on the .hove improved method is developed for three 



different modeU; a) constant properties, b) triable properties, c) thermal effects 
(energy equation) with variable properties. 

Arbitrarily varying seal profiles in both axial and circumferential directions are 
considered. An example case of an elliptical seal with trying degrees of axial cnr- 
nature is analysed In detail. A case study based on predicted clearances (6 axial 
planes with 68 clearances/plane) of an interstage seal of the SSME-ATD-HPOTP is 
presented. Dynamic coefficients based on external speciffed load are Introduced for 
seals, for the first time, to analyse seals that support a pre-load. 

The other objective of this work is to study the effect of large rotor displace- 
mmits of SSM&ATD-HPOTP on the dynamics of the annular seal and the resulting 
transient motion. Currently , the linear model of the a nn u lar seal employed at NASA 
Marshall Space Flight Center (MSFC) to estimate the seal forces during a transient 
motion of the turbopump uses a set of 6 dynamic coefficients computed at sero (s = 0) 
eccentricity. This model, while valid for a email motion of the rotor about the cen- 
tered position, may not be accurate for large off-center operation of the teal. One 
of the objectives of this study is to identify the magnitude of these deviations and 
establish limits of effectiveness of using such a model. This task is accompHshed by 
solving the balk flow model seal governing equations directly for transient seal forces 
for any given type of motion, including motion with large eccentricities. 

Based on the above study, an equivaience is established between linearized co- 
efficients based transient motion and the same motion as predicted by the original 
governing equations. An innovative method is developed to model non-linearities in 
an annular seal based on dynamic coefficients computed at various static eccentric- 
ities. This method is thorougUy tested for varions types of transient motion using 
bulk flow model results as a benclimark. 
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CHAPTER I 
INTRODUCTION 

Annular seals are used in turbomacliinery to reduce excessive leakage of the working 
fluid from high pressure side to the low pressure side. The schematic in Figure 1.1 
shows a typical seal application. The working fluid is forced to leak from a high 
pressure region to the low pressure region through the stator-rotor interface and the 
function of a seal is to reduce this leakage. Even though, their main purpose is to 
inhibit leakage, the main interest in seals from a rotordynamic point of view arise from 
the fact that the fluid forces in a seal can have a strong influence on the dynamic 
characteristics of the turbomachine, directly affecting the performance of the machine. 
Depending on a number of parameters that go into the design of a seal, these fluid 
forces may act to stabilize the rotor system, or worse, work to destabilize the system. 

Typical annular seal profiles are shown in Figure 1.2. Seals are classified based 
on the shape of the clearance profile. A seal with a constant clearance is a straight 
seal and that with a linearly varying profile is a tapered seal. A tapered seal may 
be either a convergent seal or a divergent seal depending on its slope relative to the 
direction of flow as shown in Figure 1.2. 

1.1 Literature Survey 

Extensive work has been done in the peist two decades to understand the dynamic 
behavior of seals. Nguyen (1988), provides a complete overview of the work done by 
various researchers in this area in the past twenty years. Starting with Black’s (1969) 
analysis of high-pressure seals, followed by Allaire’s (1976) eccentric seal analysis and 

Journal model is !Ibansactjons of A5AfE JovtubI of !ZWboiogy. 




straight seal convergent seal divergent seal 


Figure 1.2 Typical Seal Profiles 
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Childs’ (1985) Hits’ bulk flow analysis for tapered seeds, there has been a steady 
improvement in the modeling of annular seals and the agreement of their predicted 
behavior with experimental results. 

Of the recent work, Simon and Frene (1989) and later San Andres (1991) ex- 
tended seal analysis to include vziriable fluid properties for cryogenic applications. 
Nelson and Nguyen (1988a,1988b) are generally credited with developing the first 
finite length eccentric solution for anTnilar seals. A numerical solution is developed 
in terms of functions of Fast Fourier Transforms (FFT) which are used to compute 
the circumferential gradients of primary variables. The analysis agreed well with 
the experimental data. A improved formulation of this method is the basis for the 
current work. Yang et al. (1992) developed an ap^e^ifiiate thermohydrodynamic 
(THD) analysis for cryogenic seals. This analysis includes an approximate steady 
state solution for a centered seal. San Andres et al (1992) provided a full set of gov- 
erning equations for THD analysis based on a turbulent bulk flow model along with 
a numerical solution based on a finite difference (FDM) formulation. 

Typically, an interstage seal of the Space Shuttle Main Engine High Pressure 
Oxygen Turbopump (SSME-ATD-HPOTP) is designed either as a straight or a ta- 
pered (convergent) seal over its entire length as shown in Figure 1.2. Tests of these 
seals after a period of operation have revealed a considerable change in their clearance 
profile from their original design clearance profile. Some of these distortions are due 
to assembly and operating interferences (Scharrer and Nunez, 1989) and others are 
due to mechanical and thermal stresses acting on the housing (Scharrer and Nelson, 
1990) as the machine reaches and operates about its full power level (FPL). 

The effect of se^il distortions on rotord}maznic coefficients was first considered by 
Sharrer and Nunez (1989). They reported that a 2-D, axis3rmmetric, finite element 
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analysis which considered the internal pressure distribution, and the boundary condi- 
tions due to assembly and operating interferences produced a clearance profile which 
was wavy and different from the nominal design tapered profile. 

This distorted seal profile in the axial direction was fitted with a clearance func- 
tion in the form of a polynomial as, 

h{z) = ai -I- 02Z + osz* -I- -f a^z* (U) 

where the coeffiaents Oi, 02 , • • • etc., are coeflicients chosen to fit the distorted axial 
profile. 

They adapted the analysis of a plain seal to the case of a wavy profile seal. They 
reported a marked change in the computed rotordynamic coefficients due to a change 
in the seal profile. San Andres (1991) repeated the above study using a variable fluid 
properties model and reported similar results. 

Scharrer and Nelson (1990) treated this distortion problem using a partially 
tapered seal model. Instead of treating the distortions as a polynomial function, they 
tried to correct the predicted distortions by machining out the undesirable distortions 
at the design stage itself. The model they used to accompHsh this is a seal with a 
taper on part of length of the seal. Using this model, they conducted a parametric 
study of various performance characteristics as a function of taper length to total 
length ratio (T/L). Based on this study they recommended optimum ratio of T/L for 
best performance of these distorted seals from a rotordynamic point of view. 

Iwatsubo and Yang (1987) considered the effects of elastic deformation of the 
shaft and the seal housing due to high pressure difference, typical of a high pressure 
annular seal, and obtained dynamic coefficients based on this model. They reported 
that the direct stiffness is significantly changed when the elastic deformation is in- 
cluded. Childs (1987) studied the effects of variable radii and arbitrary clearance 
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fiinction on fluid forces developed in pump impeller shrouds. 

All the work reported in the literature on distortions in seals is limited to dis- 
tortion along the length of the seal. DetaUed thermc^tic studies based on a finite 
element model of the entire turbopump have revealed that seal distortion is not lim- 
ited to axial direction and a considerable distortion occurs along the circumference 
also. 

An example predicted seal profile of an interstage seal of SSME-ATD-HPOTP 
from a thermc^elastic finite element study of the entire pump is shown in Figure 1.3. 
In this figure, the seal is stretched out 360°, and Z is the longitudinal axis of the seal. 
The seal, initiaUy designed as a convergent seal, is severely distorted both in the axial 
direction as well as in the circumferential direction. 

1.2 Dynamic Analysis of an Annular Seal 

The mam objective of this work is to develop a dynamic analysis for liquid annu- 
lar seals with arbitrary profile and analyze a general distorted interstage seal of the 

SSME-ATD-HPOTP. The essentials of dynamic analysis of an annular seal are ex- 
plained below. 

Th. main objactiva of dynamic analysis of an annular saal is to astimata tha 
coaffiaants of tha linearized force-mation model of a saal shown in Eq 1 .2, for a email 
motion of tha rotor. Tha modal shown in Eq. 1.2 is for a two dagraa of fraadom 
(2-DOF) vibration modal. Thera are more complex models available that include 
additional dagraas of fraadom, but tha 2-DOF modal in Eq. 1.2 is tha most widely 
used one in saal literature to correlate theoretical and axperimantal data. 
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In tills equation, {Sx, 6 y) are the displacements, are the velocities and 

{Sx, 8 y) are the accelerations in the X and Y directions respectively of the center of 
the rotor, relative to a static operating point (x,y). The fluid force terms AF, and 
AF^ are the incremental or perturbed fluid forces for a small motion of the rotor 
shaft about (x,y). These force components, in general, vary as a function of rotor 

displacement, translational velocity and acceleration and are linear only for smsdl 
orbital motion of the rotor. 

In this model, F,*, k^, ky,^ Kyy are the linearized stiffness coefficients, 

C.V, Cyy are the linearized damping coefficients and M„, m^, rriy^ and Myy 

are the linearized added mass or inertia coefficients at the static operating point or 
eccentricity (x,y). 

In the linearized model, the terms [K^Jx] and [KyySy] account for the incremen- 
tal fluid reaction forces of the seal due to a small displacement of the rotor {Sx, 6y). 
The term \k^6y] is the cross coupled force in the X direction due to a displacement 
By in the Y direction. Similarly, [ky^Sx] is the cross coupled force in the Y direction 
due to a displacement Sx in the X direction. These cross coupled forces arise out of 
circumferential velocity and are a source of instability in rotor systems. The terms 
[C^Jx] and \Cyy 8 y] represent the incremental damping forces due to a small velocity 
change (8x,8y). Damping forces tend to have the opposite effect to that of cross 
coupled forces and their net effect is to add to the stabUity of the system. The terms 
[M^ 3 . 8 x] and \Myy 8 y] are the incremental fluid inertia forces due to a small change 
in acceleration [ 8 x, 8 y). For a concentric seal, K„ = Kyy, k^ = ky, etc., reducing 
the number of coefficients from twelve to six. TypicaUy, for an annular seal, the 



8 


unportant coefficients are direct stiffness, cross coupled stiffiiess, direct damping and 

direct mass. The contributions of other terms are negligible in most cases compared 
to these terms. 

These coefficients are estimated by fitting the perturbed fluid forces AF. and 

AFr due to a smaU perturbed motion about a steady state position (x,y) to the 

linearized model in Eq. 1.2. These twelve coefficients together provide a dynamic 

model of the seal for a smaU motion of the rotor and this model may be used to 

predict the fluid forces acting on the rotor in vibration-response models and rotor 
stability analysis. 

1.3 Current Work 

Th« dynamic analysis devdoped in tiis »ork is basad on a mathod originaUy proposad 
by Nalson and Ngnyan (1988a, 1988b). Thay ara araditad (Childs, 1993) noth davd- 
oping tha hrsl finila langth accantric solntion for annular saaU.Tha original analysis 
showed good agreement with experimental results. 

The original method proposed a method in which the governing nonlinear PDEs 
modeling the turbulent bulk flow in an annular seal are reduced to a set of ordinary 
differential equations by using an approximation scheme for computing the circum- 
ferential gradients of the primary variables. With this assumption, the order of the 
problem is reduced by one, i.e., from a 2-D to a 1-D problem and essentiaUy the 
problem is reduced to solving a set of ordinary differential equations for which the- 
ory IS wefl developed. This reduction in computational complexity of 

-magmtudc is the main advantage of this method compared to a 2-D finite difference 
method (FDM) or finite element method (FEM) formulation of the same problem. 

In the current work, the zeroth order and first order solutions of this original 
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analysis are unproved to make the overall solution more accurate and computationally 
more efficient and eliminate some of the reported problems with the original method. 

The original method used an innovative approximation scheme based on Fast 
Fourttr Thius/orms (FFT) to compute the drcumferential gradients of the seroth 
order variables of the eccentric solution. The numher of trigonometric functions 
included in the approximation is equal to the number of drcumferential grid points. 

as been reported (San Andres, 1991) that this method requires a considerable 
number of trigonometric fimctions for accurate solution at high eccentridties. Nguyen 
(1988) also reported problems with convergence at higher eccentricities in some cases, 
possibly due to the truncation error introduced by induding only a finite number of 
functions. In addition, another disadvantage with trigonometric functions used in the 
original method is the very CPU intensive nature of thdr computation. 

A simpler scheme based on cubic splines is found to be computationally more 
effident since it does away with trigonometric functions and it also does not require as 
many drcumferential grid points as the original method for a given accuracy tolerance. 
The increase in accuracy with cubic spline based approximation scheme ako reflects 
in better convergence properties at higher eccentridties compared to the 

method. This fact is verified by the successful analysis of cases with the new approach 
where the original method had failed. 

This flrst order solution of the original analysis is modified by induding a more 
exact solution that takes into account the variation of perturbed variables along the 
circumference. This improved analysis show better agreement with experimental 
results than earlier analysis, particularly at higher eccentridties. The new analysis 
developed treats these variables as general continuous functions and a completely 
new set of equations for dynamic analysis are derived based on this more general 
model. The original method was developed for Moody’s friction model. In the current 
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work, a unified solution procedure that is valid for both Moody’s and Hirs’ models is 
developed. 

Nguyen (1988) developed the original method for liquid seals with constant prop- 
erties and gas seals. Since the main interest of the current work is analysis of liquid 
seals for cryogenic turbopumps, the improved method will be extended to include 
variable fluid properties and thermal effects. 

Dynamic analysis based on the above improved method is developed for three 
different models. 

1. Constant fluid properties. 

2. Variable fluid properties 

3. Thermal effects (energy equation) with variable fluid properties. 

Arbitrarily varying seal profiles in both axial and circumferential directions are 
considered. The arbitrary seal profile may be either due to distortion as discussed 
earlier, or by design itself to enhance some optimum performance characteristics of 
the seal. An example case of an arbitrary profile, an eUiptical seal with varying 
degrees of axial curvature, is analyzed in detaU. An example film thickness analysis 
for this elliptical seal is presented. 

A case study based on predicted clearances (6 axial planes, 68 clearances /plane) 
of an interstage seal of the SSMI^ATD-HPOTP is presented. This predicted profile 
is obtained from a thermo-elastic finite element model of the entire turbopump. The 

results of distorted seal analysis are compared with those of a similar seal with average 
inlet and exit clearances. 

TypicaUy, seal coeflicients are computed in a minimum film thickness coordi- 
nate system as a function of eccentricity ratio and then transformed into the user 
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defined coordinate system for actual application. Such a procedure is not valid for 
an arbitrary profile seal with a non-uniform cross section. This important feature of 
directional dependence of dynamic coefficients for arbitrary profile seab is illustrated 
with reference to an elliptical seal. A method for computing these coefficients directly 
m the rotor coordinate system is presented. Dynamic coefficients based on external 

specified load are introduced for seals for the first time to analyze seals that support 
a pre-load. 

A number of cases from literature, both experimental and theoretical, arc studied 
with reference to the analysis developed. In particular, results of current work are 
compared with the theoretical work of the following: Nelson and Nguyen (1988a,1988b), 
Childs and Kim (1985), Childs and Lindsey (1993), Sharrer and Nelson (1990), Schar- 
rer and Nunez (1989), and San Andres (1991,1992). In addition, theoretical predic- 
tions from current work are compared with a number of experimental results. 

The other objective of this work is to study the effect of large rotor displacements 
of SSME-ATD-HPOTP on the dynamics of an annular seal and the resulting transient 
motion. Currently, the linear model of the annular seal employed at NASA Marshall 
Space Flight Center (MSFC) to estimate the seal forces during a transient motion of 
the turbopump rotor uses a set of 6 dynamic coefficients computed at zero (e = 0) 
eccentricity. This model, while valid for a small motion of the rotor about the centered 
position, may not be accurate for large off-center operation of the seal. One of the 
objectives of this study is to identify the magnitude of these deviations and estabUsh 
limits of effectiveness of using such a model. This task is accomplished by solving 
the hvlk flow model seal governing equations directly for transient seal forces for any 
given type of motion, including motion with large eccentricities. 

This approach of solving governing equations directly for transient seal forces 
while being the most accurate, may not be practical for routine rotordynamic sim- 
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Tilations. In fact, this is the primary reason for developing and studying approxi- 
mate linear models and their widespread use in vibration-response and rotordynamic 
simulation studies. As an alternative, an innovative method is developed to model 
non-lineanties m an annular seal based on dynamic coefficients computed at various 
static eccentricities in the seal clearance. This method, thoroughly tested for various 
types of transient motion, provides an accurate and computationaUy efficient means 
to model the effects of eccentric seal operation on the dynamics of the rotor system. 
The results from this new method compare well with bulk flow model results. 

Typically, the dynamic coefficients are computed from the first order solution 
which is directly dependent on the zeroth order solution. Even though the zeroth or- 
der equations and associated boundary conditions are essentially the same for different 
analyses based on bulk flow model, for example Childs (1985), Nelson and Nguyen 
(1988a,b), San Andres (1991), and current analysis, there appears to be variations in 
how zeroth order and first order equations are formulated and solved. In the pub- 
Hshed literature on seals, it is assumed that the dynamic coefficients extracted from 
the first order solution automatically approximates accurately the dynamic behavior 
of the original governing equations for a small motion of the rotor. Two possible 
sources that may be cited for a discrepancy between these two approaches are, a) 
maccurate formulation of the problem, b) error in implementation. In the present 
work, based on the transient analysis developed with original governing equations 
(no first order solution involved), an equivalence will be estabhshed, for the first time 
for seals, between the linearized coefficients based seal forces i.e., computing seal 
forces usmg coefficients in the linearized force-motion model of Eq. 1.2 versus the 
same forces as predicted by the original governing equations. If such an equivalence 
can be estabhshed, it proves that the dynamic coefficients being extracted from the 
dynamic analysis are indeed the correct coefficients which in turn validate the zeroth 
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and first order solutions. In other words, it is a check case for the entire analysis. 

1.4 Original Contributions 

The original contributions of the current research lire summarized below. 

1. Develop a dynamic analysis for arbitrary profile liquid annular seals based on an 
approach first proposed by Nelson and Nguyen (1988a, 1988b). The following 
modifications are incorporated into this analysis. 

(a) Improved zeroth order solution. 

(b) Improved first order solution. 

2. Dynamic analysis for eccentric seals is developed for three different models based 
on the above method. 

(a) Constant fluid properties. 

(b) Variable fluid properties. 

(c) Thermal effects (energy equation) with variable fluid properties (concentric 
case). 

3. A unified solution procedure is presented for the following two friction models. 

(a) Moody’s Model 

(b) Hits’ Model 

4. Dynamic coefiicients for seals based on external load specification. 

5. Application of the new method to study the static and dynamic characteristics 
of an arbitrary profile seal, e.g., an elfiptical seal with a varying axial curvature. 
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6. Study of directional dependence of dynamic coefficients for arbitrary profile 
seals. 

7. Dynamic analysis of a general distorted interstage seal of SSME-ATD-HPOTP 
turbopump. 

8. Tyansient analysis with an annular seal for large eccentric motion of the rotor. 

(®) Tr&usient analysis with bulk flow governing equations. 

(b) Comparison of bulk flow model simulations with linear model (dynamic 
coefficients computed at concentric position) results. 

(c) Study of equivalence between rotordynamic coefficients based transient 

motion and the same motion as predicted by the original governing equa- 
tions. 

(d) A new method to model non-linearities in an annular seal for transient 
analysis. 

(e) Thorough testing of the new method for various types of transient motion 
and comparison with bulk flow model. 



15 


CHAPTER II 

CONSTANT PROPERTIES MODEL 
2.1 Bulk Flow Governing Equations 


The bulk flow governing equations of turbulent fluid flow in an annular seal have 

been derived using several approaches. The following analysis is based on the work 
of Nelson (1984). 

The primary variables of the bulk flow are the axial velocity 11(2,5), circumfer- 
ential velocity v{z,q) and pressure p{z,q). The variation of these primary variables 
across the thin film is neglected. The axial and circumferential coordinates are z and 
q respectively. The radius of the rotor is R and the rotor angular velocity is u; rad/s. 
The differential fluid volumes used to derive the bulk flow governing equations 

are shown m Figures 2.1-2.2. Mass conservation of fluid in the differential volume of 
Figure 2.1 yields, 

p{v+.^ig)(h+-g-dq)dz - phvdz + p^„+^dz)^h+~dz)dq - puhdg+p^dzdg = 0 


( 2 , 1 ) 


or. 


dh du dh dv dh 


dz ' ''dq ^ ^ dt = ^ (2.2) 

Conservation of momentum in the axial direction for the fluid in the differential 
volume of Figure 2.2 may be expressed as. 


(phdzdq) - ( r„-r„)dzd5 + phdq - {jp^^dz){h-\-^dz)dq -|- P^dzdq (2.3) 

where The terms are rotor and stator surface shear stresses in the axial di- 

rection. ^ is the total or material derivative of the axial velocity u and is defined 



16 


as. 


,du dll du. 

Dt ~ 


Simplifying the above expression yields, 


— ^ 9u 1 

dz ~ + + + + ’•'■) (2.6) 

Similarly, conservation of momentum in the circumferential direction is 
. Dv 


,du du 


dz ' dq‘ 


(2.4) 


is expressed as, 


{phdzdq) ( r„ Tr,)dzdq + {p ^dq){h ^dq)dz - phdz (2.6) 

where r,„r., are the rotor and stator surface shear stresses and ^ is the total deriva- 
tive of circumferential velocity t; defined as, 


fdv dv dv. 

m - + 


or, 


, ov av dv 1 

dq dt'^'^Tz'^ '’di^ + 

The shear stresses at the wall based on Moody friction factor are given by, 

-K« + TV,) = p/,^v/u2-t-v2 -h Pfr^^~^\/V-^^{V-WY 
-(r„-fr„) = p/,^v/n2-fv2 -f pfr'^yju'i ^ [v - wY 

where, 

P density of fluid 

w rotor rpm in rad/s 

^ wii, rotor surface velocity 

ft stator friction factor 

fr rotor friction factor 


(2.7) 


( 2 . 8 ) 


(2.9) 

( 2 . 10 ) 
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(tt + ^dz){h + ^dz) dq 


Figure 2.1 DifFerential Volume for Deriving the Continuity Equation 



Figure 2.2 Differential Volume for Deriving the Momentum Equations 
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Using the transformation for the circumferential coordinate q = R/3, m the 
Eqs. (2.2, 2.5, 2.8) yield the following bulk flow continuity, axial momentum and cir- 
cumferential momentum equations for an incompressible fluid. 

Continuity: 

1 a(hv) ah 

dz + R-af + « = " ( 2 - 11 ) 

Axial Momentum: 

, « du du^ 

pdz ^dt Rd/3 

+ + V* -I- fr-yju^ + (v - wY (2.12) 

Circumferential Momentum: 

V dv dv, 

pRd/3 ^dt R~^ 

+ f.^y/u^ + v^ -f- ^2.13) 

2.2 Friction Factors 

Two friction models extensively used in seal analysis are the Moody’s model and 
the Hirs’ model. These two models differ in the way the roughness of the surface, 
both stator and rotor, are modeled. Of these two models, use of Moody’s model 
is more prevalent because of a more reaUstic friction factor which is dependent on 
local Reynolds number, film thickness and surface roughness compared to the Hirs’ 
model where the coefiicients are for a fixed clearance and an average Reynolds num- 
ber. However, considerable experimental and theoretical data exists for Hirs’ model 
based analysis, for example ChUds (1985), Sharrer and Nelson (1990) etc., making it 
attractive for comparative studies for any new analysis such as the current work. 
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Moody’s Model: 


f. 


fr 


R, 

Rr 


0.0055 

4 

0.0055 

4 


{ 1.0 + 
( 1.0 + 


2ph 

-j- 


( 10 ‘| + 



10 « _)./>} 

10«±)V3} 


where, 

stator pocket 
Cr rotor pocket 

h film thickness 

c* nom in al radial clearance 

R§ stator Reynold's number 

Rt rotor Rejm old’s number 

2^ stator relative roughness 

rotor relative roughness 


Hirs ’ Model: 


f. - 

fr = 

where B„ m„ n. and 771, are Hire’ constants for stator and rotor respectively. 


(2.14) 

(2.15) 

(2.16) 
(2.17) 


(2.18) 

(2.19) 
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2.3 Film Thickness 

Tie expreetion for film thickness k{z,0) as a function of eccentricity is derived in 
the rotor (fixed) coordinate system, instead of a mimmum film thictma, coordinate 
system notmaUy used in straight or tapered seal analysis. In a typical analysis with 
these seals, the minimum film thickness coordinate system is usually aligned with 
the X-axis of the rotor coordinate system and for eccentric operation, the dynamic 
coeffiaents ate computed as a function of eccentricity along this axis. Although the 
dynamic coefficients can be rotationaUy transformed once they are evaluated at a 
given equffibrium position, the minimum film thickness system cannot be used with 
a seal that has a circumferentially varying clearance (which destroys axisymmetry). 
This same restriction applies to, for instance, modeling pressure dam bearings versus 
plain journal bearings. As a result of this asymmetry, the dynamic coefficients have 
to be computed directly in the rotor coordinate system for non-uniform profiles or if 
a minimum film thickness system is used, the orientation needs to be specified. This 
important feature of the directional dependence of dynamic coefficients for arbitrary 
profile seals will be further discussed with reference to an elliptical seal in Chapter V. 

The seal geometry is, in general, defined by its clearance function c(r,/3). The 
clearance function of a seal defines the fluid film thickness when the rotor is at the 
centered or concentric position with respect to the seal. A constant c specifies a 
straight seal, a linear function in r defines a tapered seal and so on. For the purpose 
of this study, any profile other than a straight or a tapered profile is considered as 
an arbitrary profile. The film thickness, which varies with eccentricity, is derived as 
a function of c(z,0) and the eccent.ricity vector (e,*). The angle if,, defined as the 
eccentnaty angle, is the angle made by the eccentricity vector with respect to the 
fixed X-axis. The magnitude of the eccentricity vector is given by the eccentricity. 
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Figure 2.3 Diagram for Deriving General Seal Clearance Expression 

e. The expression for the film thickness is given below with reference to Figure 2.3. 
In this figure, e*, Cy are the offsets of the center of the rotor O' with respect to the 

center of the seal denoted by O. In this figure, <f> is the eccentricity angle and is 
the angular coordinate. 

By the law of cosines. 


or. 


(i? + h)2 = e* + {R + cy -2e{R + c)co8{fi-<f>) 


( 2 . 20 ) 


= \/e* + (i2 + c)2 - 2e(R + c)cos(fi - </,) - R 


( 2 . 21 ) 


ecos{fi — <i>) — egCos/3 + tysinfi 


substituting. 


e, 


ecoa<f> 


( 2 . 22 ) 

(2.23) 



22 


= esin<}> 


(2.24) 


Eq. (2.21) may be rewritten as, 

= \/(^ + c)2 - (e.sin/3 - tycosfiy - (e,coa/? + e^atn/?) - iE (2.25) 

and its gradients in /? sind z directions are, 

^ _ (-R + c)^ — [e^sinfi — eyCoj/?)(egC03y3 -f Cysi nff) 

cY — — tyCosfiY 

+ {e^sinfi — eyCos/3) (2.26) 


5/i 


(^ + c)|f 


^(/Z + c)2 - (e,stn/3 - e^cos^Y 

Besides specifying the film thickness in a fixed coordinate system, the above 
expression for film thickness has the following advantages over the more commonly 
used approximate form. 


— c(z) - e^cosfi — CySinY (2.28) 

1. It models the curvature of the film thickness accurately. This is important, 
particularly, when analyzing a severely distorted seal or an arbitrary profile seal 

with a clearance function varying in the circumferential direction, such as an 
elliptical seal. 

2. It specifies the film thickness in a fixed coordinate system which is essential 

for analyzing non-uniform profile seals as explained earlier, (examples to be 
discussed later). 

3. The general expression in Eq. (2.25) is more accurate, mathematically, partic- 
ularly at higher eccentricities than the approximate form in Eq. (2.28). 
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2 A Solution Procedure 

The general steps involved in the solution procedure are outlined in the flow chart 
(Nguyen, 1988) shown in Figure 2.4, These various steps are summarized below. 

1- Derive the bulk flow governing equations. 

2. Perform perturbations on the original governing equations to yield zeroth order 
and first order governing equations. 

3. Form the appropriate boundary conditions at the inlet and the exit. 

4. Solve the set of zeroth order equations subject to the boundary conditions at 

inlet and exit to obtain the zeroth order (steady state) solution of the primary 
variables, uo, vq, po- 

(a) Compute leakage. 

(b) Compute steady state reactive seal forces. 

(c) Compute frictional torque. 

5. Perturb the zeroth order boundary conditions to obtain the first order boundary 
conditions. 

6. Assume a harmonic solution form and use a separation of variables procedure 
to reduce the first order equations to a set of ordinary differential equations. 

7. Solve for the first order variables, u^, ui, p„ subject to the first order boundary 
conditions. 

8. Extract dynamic coefficients from the first order pressure field. 



Governing Equations 



Figure 2.4 Flow Chart of Solutiou Precodure, Nguy«. (1988) 
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2,5 Perturbation Analysis 

In this section, the zeroth order and first order equations are derived using a per- 
turbation analysis. The original bulk flow governing equations are perturbed about 
their steady state values to yield the zeroth order and first order governing equations. 
The zeroth order equations are also known as the steady state equations and they 
may also be obtained from the original governing equations by dropping the time 
dependent terms. The perturbed or first order equations of the governing equations 

Eqs. (2.11-2.13) are derived for a small motion of the rotor about a steady state 
eccentric position. 

The assumed form for the dependent variahles and film thickness for perturbation 
analysis are given €is, 


u[z,0,t) = Uo{z,f3) -f- e«i(z,/3,<) (2.29) 

v(z,/3,t) = vo{z,P) -I- eui(z,/3,<) (2.30) 

p{z,/3,t) = po{z,0) -f epi{z,^,t) (2.31) 

h{z,/3,t) = ho{z,0) -f ehi{z,/3,t) (2.32) 


where c is a smaU perturbation and uo, vo, Po, ho are the zeroth order variables and 
«i, vi, Pi, hr are the corresponding first order variables. Substitution of these expres- 
sions into Eqs.(2.11-2.13) and neglecting second and higher order terms yields sets 
of zeroth order and first order equations of the form, 


(zeroth order equations) + e(first order equations) = 0 
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2.5.1 Zeroth Order equations 


Tht zerotl order equations are essentially steady state equations and they may also 

be obtained by canceling out the time dependent terms in the original governing 

equations. The subscript 0 in the following seroth order equations refer to zeroth 
order variables. 


h^dpo 
Po dz 


d{hoUo) 1 5(^Vo) 

dz R dp ^ 


ho{ 


Vo duo 

rW 


+ 


, duo. 


fio ^ -T /.0-^^Uo + (no — tn)* 


(2.33) 


(2.34) 


^Po _ , jVodvo dvo. 

PoR dp - ^°^Rjp + 

+ /.o-j\/«o + Wo + fro ~ ^ y/ul 4- (tfo - w)^ (2.35) 


where the friction factors U and U are the friction factors. 

The set of equations in Eqs. (2.33-2.35) are the zeroth order equations for a 
constant fluid properties model or an incompressible fluid. 


2.5*2 First Order Equations 

The first order equations are given by, 

Continuity: 

^^0 _ ^^0 Vo dho 

dz ~ R~^ 

,du 1 


5uo 


+ 


ho dvo 

'rTp 


+ 


^0 

dz 


■Ui + 


1 dho 
R dp 


Vi = 
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Axial Momentum: 


L ^0 ho dpo 

houo^ + —p- 
Oz po dz 


+ 


hovo duo ^ duo 


+ A^ui + A^vi = Ahhi (2.37) 


Circumferential Momentum: 


h . 

«o«o-^ + 


hoVo dvp ^ ho dpo ^ dvo 

E d/3 


Bnhi (2.38) 


where the coefficients A^, A„, Ah, B„, B„ and Bh are functions of steady state vari- 
ables «o, Vo, Po and their axial and circumferential gradients and friction factors and 
their derivative terms. These coefficients for constant properties model are given in 
Appendix A. 

It can be seen that the first order continuity and momentum equations do not 
change between the Moody’s and Hirs’ friction models and the friction factor model 
only affects the definition of the coefficients >1, • . . etc.. These coefficients are 
derived m such a form such that the solution procedure is valid for any general friction 
model. Specific analyses for two particular models, a) Moody’s model b) Hirs’ model 
are developed based on this general format. 


2,5.3 Linearization of Ftiction Factors 

The friction factors /, and /,, for constant fluid properties, are implicit functions of u, 
V and h. The perturbation in the friction factor is obtained by a linearization process 
using Taylor’s series expansion about the operating point. The foUowing example 
analysis illustrates the steps involved in the linearization of friction factor /,. Using 
Taylor’s series expansion of /. about the steady state variables, uo, vo and po, 

ft{u,V,h) = /,o|(uo,«o.fco) + ^l(•oJ>o)(« - «o) 

+ §7l(«.«.)(f--yo) + ^\(^.vo){h - ho) (2.39) 
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or. 




(2.40) 


The expressions, are derived for both Moody and Hirs friction factor 

models and are given in the Appendix D. 


where, 




_ Uo 



du 

^•0 2 1 2 

(2.41) 


d_u 

_ Vo 



dv 

3*0 2 1 2 

«o + 

(2.42) 


df. 




dh 

ho 

(2.43) 


0.0055 

X 10® .t>. 1 


9to = 

I2ie„ + 

(2.44) 

hto = 

0.0055 

12 

(10^^ + 10«-^)'/» 
ho R^' 

(2.45) 

RfO = 

2/>Q^Q 

Ao + Vo 

(2.46) 


For the case of a fluid with variable properties, there wiU be two additional 


Hi Hl 

ap’Bii’ 


terms. 


2.6 Zeroth Order Boundary Conditions 

The boundary conditions for the zeroth order or steady state equations are illustrated 
in Figure 2.5-2.6. 

The fluid flow in an annular seal occurs from the high pressure side (inlet) to the 
low pressure side (exit) as shown in Figure 2.5. Just prior to the inlet {z = 0), the 
fluid has zero axial velocity and the fluid pressure is given by the supply pressure or 
reservoir pressure p. . At the entrance of the seal a swirl is induced in the fluid by the 
eye of the impeller creating the tangential velocity of the fluid as shown in Figure 2.4. 
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GcneraUy, the magnitude of this tangential velocity is estimated as a percentage of 
the rotor surface speed and is specified by the prc-5u«W ratio as, vo = psr{u;R), where 

R is the radius of the rotor, u is the angular velocity in rad/s and psr is the pre-swirl 
ratio. 

At the inlet, as the fluid enters the seal there is a loss in pressure with a corre- 
sponding increase in the acceleration of axial velocity, u. The relationship between 
these two variables is given by the Bernoulli’s equation as, 

Pi - Poi(0,/3) = -po«w(0,^)(l -l-fi) (2.47) 

where, is the inlet loss coefiicient and poi, uoi refer to the pressure and axial velocity 
right after inlet as shown in Figure 2.6. 

The inlet or entrance loss factor, in general, is a function of geometry at the 
entrance as well as local Reynold’s number (Nguyen, 1988). In the present work, 
as is the norm in seal literature, a constant value is assumed for this coefficient. 
In practice, one of the methods used to estimate the pre-swirl ratio and the inlet 
loss coeffiaent is by matching the theoretical flow rate with experimentally measured 
data. At present, there is no rehable way of predicting these two parameters and 
the empirical procedure used above typically gives rise to different sets of input data 
depending upon the seal analyst’s objectives. 

Right after the exit {z = L), the pressure is given by the low pressure, p,. At 
the exit, a similar relation as given in Eq. (2.47) is used to relate the variables. 

PO2(0,/?) - Pe = 2^ottM(0,^)(l -^,) (2.48) 

where, e, is the exit pressure recovery coefficient. Typical values for a worn and new 

seal are 0.7 and 0.85. In this analysis, the value of = 1.0 is used, i.e, there is 
complete recovery of pressure. 
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The following is a summary of the boundary conditions for the zeroth order 
equations. 

At the inlet: 

axial velocity, tt©: 

prior to inlet: 

tto(0,/3) = 0 (2.49) 

right after inlet; 

uo(0,/3) = uoi(0,/?) (2.50) 

circumferential velocity, v©: 

vo(0,/?) = psrXLoR (2.51) 

pressure, p©: 
prior to inlet: 

Po(0,/3) = Pi (2 52) 

right after inlet: 

Po(0,;5) = poi(0,^) (2.53) 

The pressures pi, poi{O,0) and axial velocity u©i(0,^) at the inlet are related by, 

Pi - Poi{Q,/3) = ^Po«oi(0,y3)(l + ^i) (2.54) 

At the exit, the exit pressure recovery coefficient is assumed to be 1, i.e., 

Po2{O,0) - pe = ( 2 . 55 ) 

or. 


P02{0,I3) = Pe 


( 2 . 56 ) 
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Let Ap be the total pressure drop across the seal. 

Ap = Pi -p, (2.57) 

Eq. (2.54) may be rewritten as, 

Poi(0,/?) = Ap +p, - ^po<(O,;0)(l + e.) (2.58) 

For the case of an incompressible fluid, the absolute pressure is not important 
and therefore the exit pressure may be set to zero, i.e.. 


Pe = 0 

(2.59) 

Eq. (2.58) may be rewritten as, 


Poi(0,/3) = Ap - ^po<(0,/?)(l + ej 

(2.60) 

or 



(2.61) 

At the outset, poi(0,/?) is unknown and must be solved iteratively by requiring 

that the pressure distribution at the seal exit satisfies the following 

condition. 

• 

K> 

II 

o 

(2.62) 


subject to the constraints of Eqs. (2.51) and (2.54). 
2.7 Solution Procedure for Zeroth Order Equations 


The solution for zeroth order equations involves the direct integration of three coupled 
nonlinear partial differential eqaations subject to the bonndary conditions given in the 
previous section. The current analysis uses a leroth order solution procedure different 
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from the original method. The following is a brief comparison between Nelson and 
Nguyen s original approach and the one used in current work. 

2.7.1 Comparison with Nelson and Nguyen Approach 

The three steady state equations are arranged in the following fashion and integrated 
from inlet to the exit. 

^ = J^v(uo,vo,po,^,^,^) (2.63) 

, ^ . ^p(^o^^o,Po, 

The functions F„, F„, Fp, for a constant properties model, are given in Appendix 
A. 

The original analysis of Nguyen (1988) proposed a method by which the coupled 
partial differential equations are reduced to coupled ordinary differential equations by 
approximating the circumferential gradients of the variables vo and po as shown in 
Eq. (2.63). At each axial step in the iterative procedure, the gradients with respect 
to P are computed based on the values of the variables at the previous step. An 
approximation scheme based on Fast Fourier Transforms (FFT) was used for this 
purpose. Assuming that velocity and pressure distributions are known at an axial 
step, the gradients of these variables, may be computed by specifying 

each of the known values as a finite-length complex Fourier series as given below. 

“o(i^j/?) = Real{uo + 2^ tt„(z)e"*^} (2.64) 

Vo{z,P) = Real{vo -f 2 j Vn(z)e*"^} (2.65) 

Po{z,P) = Real{po + 2 j p„(z)e‘"^} (2.66) 

where N is equal to one half the number of circumferential divisions. The above 
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expressions may then be used to compute the gradients as, 

duti{z,3) fN-i 

^ = Real{2j^^ mu„(z)e”^} (2.67) 

dvo(z,/3) fN-l 

= Real{2j^^ mv„(z)c*"^} (2.68) 

dpo(z,3) rAT-l 

inp^zy”^} (2.69) 

This methods suffers from the following drawbacks. 

1. It requires too many functions for accurate solution at high eccentricites (San 
Andres, 1991). 

2. Computation of trigonometric functions is a very CPU intensive procedure. 

3. Convergence probems (Nguyen, 1988), possibly due to truncation error intro- 
duced by not including enough functions in the approximation scheme. 

4. No rehable way to decide on the number of functions to be used. 

In the present analysis, a simpler method based on cubic splines is implemented. 
This method is more accurate as no truncation error is involved as in the FFT method. 
Also, convergence at higher eccentricities is achieved with relatively fewer iterations 
than the FFT method. It is also computationally more efficient as it does not involve 
the computation of CPU intensive trigonometric functions. Also, the number of 

circumferential grid points may be varied, upto a limit, with out affecting the accuracy 
of the solution. 

2.7.2 Iterative Solution for Zeroth Order Equations 

TypicaUy, an iterative procedure is used to solve for the pressure distribution. Fig- 
ures 2.7-2.8 illustrate typical subdivisions in the axial and circumferential directions 
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for spline interpolation and numerical integration. The circumference is divided into 
segments of equal length both in the axial and the circumferential directions. 

Let m be the number of axial grid points, while n is the number of circumferential 

grid pomts. Let the superscript j indicate a function evaluated at /? = fij. The coupled 

partial differential equations are reduced to a set of coupled ODE’s, by moving all 

^-dependent terms to the right-hand side as shown in Eq. (2.63). At any given axial 

plane, cubic splines are used to lit and to compute the gradients 

&vn (^) 1 &pfj (J) V * 1 ft 

^0 B 0 where j = 1,2, • • • ,n the number of circumferential grid points, at 

that axial plsme. 

The set of differential equations in Eq. (2.63) are solved using separation o/vari- 
ables type assumption that at any fixed axial location, the gradients of the dependent 
vanables uq, vo, po are obtained by spline fitting these variables at that axial location 
as a function of the circumferential coordinate (3. 

Let e represent a cubic spline function operator which when applied to a set of 
function values yields a continuous piecewise cubic spline approximation in /3 of that 
fimction. In other words, given a set of grid tmiues for ,8 = ft, ft, . . . , nt any nrinl 
location z = Zk, 


Mzk,/3) = J = 1, 2, ..., n (2.70) 

vo{zk,fi) = 0(/?,t;{^)), j = 2, ..., n (2.71) 

= Q(/3,Po^), j = l, 2, •••, n (2.72) 

The circumferential gradients, g, & are evaluated as, 

_ a0(/?,4")), 

dfi ~ (2-73) 

d/3 d/3 


(2.74) 




Fignre 2.7 Circumferential Grid 



Figure 2.8 Grid for Numerical Integration 
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d/3 d0 (2-75) 

The set of ODE’s are solved subject to the boundary conditions 
At the inlet, (z = 0), 

circumferential velocity, 

= psr X (wiZ), j = 1, 2, •••, n (2.76) 

pressure, 


+ = 1, 2, ..., (2.77) 


At the exit, {z = z„,), 


Po'^(0,/?i) = 0, j = l, 2, ...,n (2.78) 

In the parlance of numerical analysis, the above problem is classified as a “ two 
pomt boundary value problem”, since the known boundary conditions exist at both 
ends of the boundary. In this study, this problem is solved using a multi-dimensional 
Newton-Raphson method known as “shooting method”. Nguyen (1988) reported 
using a similar method based on numerically computed gradients. 

The problem may be specified in the following terms. It requires to find the n 
unknown inlet pressures, p[ ^(zj), p^^^(zi), •••, p^**^(zi) subject to the condition that 
the n outlet pressures pS^^(z,„), pfW), • * *, P^e\zr.) are forced to be zero. 
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This may expressed sls, 

Po ^(2l),Po 0 

Po ^{■*TTnPo ^(■2l)>Po ^(^1)5 ■ ' • »Po”^(^l)} 0 

= < [ (2.79) 

Po ^{^>Po ^(•2 i)jPo '(2^i)> ■ ' • jPoTH^i)} 0 

It may be noted that even though both and are unknown at the inlet, 
there are only n unknowns at Zi due to the relationship given in Eq. (2.77). 

The two point boundary value problem is now reduced to finding the n roots 
of the above equation. The numerical integration of the preceding differential equa- 
tions may be considered as only a means of evaluating the functions Po\z„t) for any 
given set of guesses for the roots Pc/^(ri), po ^(^i)j ’ • As mentioned earlier, 

shooting method is an application of a multidimensional Newton-Raphson method 
for iterative search of roots Pq ^(ri), Pq ^(.si), •••, Po”^(i^i) combined with a numerical 
I^^^csr^tion based on evaluation of the exit pressures Po^'^(zm) 

Dropping the subscript 0, let pl^^(zi) represent the k-th guess for the set of inlet 
pressures Po^^(zi) that satisfy the boundary condition of Po’^Zm) = 0 at the exit. Let 
the next guess be given by 


p1+i(«i) = PkH^i) + Apfc'^(^i) 


(2.80) 
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where Pfc+i)( 2 i) is the new guess. This may be expressed mathematically 



(2.81) 



The partial derivatives are computed numerically using a finite dilFerence formula, 
^Po^(gm) ^ 1 

bo^{^m>Po ^(^l),-- *,Po'^(2l) + Ap^-'^(zi),- . . ,p^”^(zi)} 

~ Pok^jPo ^(^i)>'‘*tPoi^^(i:i),---,Po"^(zi)}] 

J ~ 2, ' • • , n; t = 1, 2, • • • , n (2.84) 


It requires one complete numerical integration of the Eq. (2.81) to compute 
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P'k\zm) in Eq. (2.82) and an additional n integrations of Eq. (2.81) to obtain the 
denvatives in Eq. (2.83). These integrations may be performed by any standard 
integrators. Nguyen (1988) used the simple Euler’s method for the above numerical 
integration. For this study, the following are used. 

1. 4-5th order Runge-Kutta-Fehlberg Method. 

2. Predictor- Corrector Method. 

3. Adams’ Method. 

2.7.3 Cubic Spline Interpolation 

Cubic sphnes are widely used in interpolation and surface fitting and a brief to intro 
duction to splines is given in this section. 

Let the dependent variable pressure p, along the circumference of the seal at a 
given axial location z = be specified by p(z„,/3,) = P(/?^.), for j = 1, 2, • • • , n 

corresponding to /? = /3j, . . . , ^^ere n is the number of circumferential grid 

points. 

In a given mterval (/?^,/3j+i), the dependent variable P{(3) is interpolated using 
a cubic polynomial function of the form given below. 

P{0) = Pi + M/?-/?i) + d,(/?-/?,f, (Si<p<ISi^, (2.85) 

The above equation may be rewritten in a more general form to facilitate the 
computation of the linear coefficients 6^, c,-, d,. These coefficients are determined 
such that the above cubic polynomial is reduced to a cubic spHne function, i.e., the 
function P(/3) and its derivative P'(/3) are continuous in (/3,,;3,^,). In the foUowing 
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equation " refers to the second derivative of the dependent variable. 


PW = 

APi + BPj^y -f CP': -f DP':^, 

(2.86) 

where, 



A 

_ i^j+i — 0 


/i 

Pj+l ~ Pi 

(2.87) 

B 

_ P- Pi 
Pi+i — Pi 

(2.88) 

C 


(2.89) 

D 


(2.90) 

The above formulation hsis the following important features. 



1. The function P(/3) is continuous at fij and i.e., at Pj in the interval 

iPj-iiPj) at in the interval 

2. The function P{f3) has a continuous derivative at and simn^r to the 
above condition. 


3. The above two conditions give rise to four constraints which are used to deter- 
mine four adjustable linear coefficients in terms of P,- P^ . , P" P" 


4. Also, four is the number of parameters required to define a cubic polynomial 
general. 

The first and second derivative of the interpolating function is given by, 

3A*-1 


m 


^ - Pj 

d/3 

3P*-1 


+ 


6 (/?i+i-^i)i^;i 


— - AP" + BP“ 


(2.91) 

(2.92) 
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By imposing the condition of continuity in first derivative at node points, the 
following equation is obtained for j = 2, 3, • • • , n - 1. 

tzp=lp; + p" + - Pj Pi - Pi-, ,, ,,, 

This is a system of n — 2 simultaneous linear equations involving n unknowns, 
2, • • • , n. Two additional conditions are required to uniquely define the 
cubic spline. These two conditions may be provided in various ways based on the end 
conditions i.e., at f3i find /3„, 

Let si{/3) and s„{/3) be two cubic polynomials that pass through the first and the 
last four data points. The two end conditions that complete the solution are defined 
by forcing these two cubics to have the same third derivative at the end points. 


^'"(A) = P"'{0i) (2.94) 

s"'{0„) = P“'i0r.) (2.95) 

The £^() may be rewritten in a more simplified form for actual computations. 


A/3,- 

1 

II 

(2.96) 

A.- 

_ Pi*, -Pi 

(^ 3+1 — 03 

(2.97) 

Af 

_ A,+i — A,- 
03+i ~ 03 

(2.98) 

Af 

_ A;«-Af 

03+^ ~ 03 

(2.99) 


and the end conditions, 


« (A) = 

6Ai’> 

(2.100) 

> {0n) = 

«aS, 

(2.101) 



43 


The quantities, A,,2A^*\ 6Aj®^ are approximations of the first, second and third 
derivatives respectively. 

-A/3i A^i 0 0 ... 0 

A^i 2 (A/?i + A/32) A/32 0 i 0 

0 A^2 2(A/?2 + A/?3) A/^s : 0 

• • 

• » 0 2(AA._J + A/3„_, 


' p; ' 



p" 


A2 ~ Ai 


^ ^ 

A3 — A2 

^'-1 


An-l “ An-2 



-A«_.a 23 


( 2 . 102 ) 


The above system of equations has the Mowing important characteristics. 

1. The matrix is diagonal. 

2. The matrix is symmetric. 

3. The matrix is nonsingular and tridiagonal. 

4. Efficient matrix reduction techniques available for solution. 


The ongmal cubic polynomial is given below along virith the coefficients in terms 

ofPi,p,+i,p;,p;;,. 

Typically, the cubic spline is written in the following form. 
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In the interval (/?,•, /3,+i), 

P{P) = Pj + fci(/3 - (3j) + cj{0 - (3jf + dj{P - fij)\ ^j<0< /3j+i (2.103) 




3P; 

P" — P“ 

^j+l 

A/3,- 


(2.104) 

(2.105) 

(2.106) 


The sets of coefficients are stored for the n — 1 intervals and the function values 
and its derivatives are computed whenever they are needed. 


P'{(3) = 6i + 2c,(/3 -/?,) + 3d,(/?-/3,f, (2.107) 

P"{0) = 2c,- + 6d,-(/?-^,-) (2.108) 

/?i < /? < /?i+i 


2.7.4 Leakage 

The mass flow rate, dQ, through the differential element h^Rdfi shown in Figure 2.9 
is given by, 

dQ = uo(0,/3)/>oho(0,/?)iZd/? (2.109) 

where, 

ao(0,/9) axial velocity at inlet 

ho(0,/?) film thickness at inlet 

R radius of the rotor 

Po density of fluid 

The above expression is integrated around the circnmference of the seal at the 
inlet to give the total meiss flow rate or leakage, Q. 
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Figure 2.9 Leakage 



uo(0,/3Vo/k)(0,^)iW/3 


( 2 . 110 ) 


2.7.5 Steady State Seal Forces 

The reactive seal forces acting on a non-vibrating rotor are obtained by integrating 

the zeroth order pressure field, po(z,/3), around the rotor and along the length of the 
seal. 

The X and V components of force acting on the differential area element Rdfidz, 
shown in Figure 2.10, are given by, 

-dF^ = po{z,fi)cosfi RdfSdz (2.111) 

-dFy = po{z,0)sin/3 Rdfidz (2.112) 

Integrating the above force expressions over the entire surface area of the rotor 
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yields, 


Figure 2.12 Frictional Torque 


-F^ 



Po{z,l3)cosfi Rdf3dz 


(2.113) 


rL 

~ ~ Jo Jo Rd0dz (2.114) 

The angle made by the resultant seal force F with the Jf-axis is defined as the 
load €mgle €md is given by, 


V- = tan ‘(-^) (2.J15) 

P = ^/FF+^ (2.116) 

The resultant seal force is abo known as the load bearing capacity of the seal. 
This resultant seal force must balance against the external load (pre-load) appHed 
by the rotor on the seal. The pre-load may vary in magnitude and direction as the 
pumps’s speed or power-level changes. 


2.7.6 Fkiction Loss 

The friction loss or horse power loss computation is illustrated in Figure 2.12. 

The frictional torque on a differential element Rd^dz due to friction at the rotor 
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2.8 First Order Boundary Conditions 

The first order boundary conditions are obtained by perturbing the zeroth order 
boundary conditions of Eqs. (2.51,2.54,2.56). 


Pi(0,/3) = -(1 + ^,)pouo(0,/3)ui(0,/?) 

(2.121) 

vi(0,y9) = 0 

(2.122) 

Pi(I,/3) = 0 

(2.123) 


2.9 Solution Procedure for First Order Equations 

The set of first order equations in Eqs. (2.36-2.38) are further reduced by employing a 
separation of variables technique for an assumed small motion of the vibrating rotor. 
The assumed form of perturbations for the dependent variables and film thickness 
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Figure 2.13 Elliptical Whirl Orbit 


are, 


u 

= Uo 

+ 

eul 

(2.124) 

V 

= Vo 

+ 

evi 

(2.125) 

p 

= Po 

+ 

epi 

(2.126) 

h 

= ho 

+ 

ehi 

(2.127) 


The rotor is assumed to execute a whirling motion with an elliptical orbit as 
shown in Figure 2.13. The figure shows the steady state operating position of the 
rotor along with the vibrating (whirling) rotor. Let the semi-major and semi-minor 
axes of this perturbation ellipse be given by (X,Y). 

The perturbation c may be considered to be a combination of two individual 
perturbations e, and e^. This assumption is the basis for eccentric seal analysis 
where the dependent variables vary in the circumferential direction as opposed to 
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the concentric case where they remain constant in the circumferential direction. The 
magnitudes of these perturbations are arbitrary since they do not form a part of the 
final solution. These perturbations are assumed to be infinitesimally small and is 
the basis for neglecting second and higher order terms in the perturbation analysis. 

These perturbations, c, and are the non-dimensionalized axes of the whirling orbit 
as given below. 


eui 

X Y 

= — + — Uly 

c. c, ' 

(2.128) 

evi 

II 

(2.129) 


X Y 

(2.130) 

ehi 

“h 

c. c* 

(2.131) 

where c. is some nominal clearance used to non-dimensionlize the perturbations. 

Let the perturbations be redefined 

in the following non-dimensionalized form. 


> 

II 

(2.132) 


Ae, = I 

c. 

(2.133) 

Substituting Eqs. (2.132-2.133) in 

Eqs. (2.128-2.131) yields, 


eui = 

AC,Ul, + ACytljy 

(2.134) 

€Vi = 

ACeVix + ACyVly 

(2.135) 

epi = 

ACxPl, + AtyPlj, 

(2.136) 

eh, = 

Ae^hig Acyhiy 

(2.137) 


Assuming that the rotor whirls about its equilibrium position in an elliptical orbit 
whose center is located at (xo,yo), then the position of the center of the vibrating 
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Figiire 2.14 Perturbation Orbit 

rotor relative to its static eccentric position is given by (Figure 2.14), 


X — Xq = Xcosutt 

y ~ Vo = Y cosu)t 


(2.138) 

(2.139) 


Let Q _ oit, where w is the angular velocity of the rotor given in rad/s. The 
perturbed film thickness expression may be rewritten as, 


h[z,(3,t) - ho{z,0) - {x-xo)coafi - (y - yo)am/3 (2.140) 


or, 

= ho{z,P) + €hi(z,/3^t) (2.141) 

Prom the above equation, the perturbation in film thickness ehi is given as. 


— —Xcosacos/3 — YainctsinP 


(2.142) 
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or 


hr 

dhi 

d/3 

dhi 

dt 


~~{Ae* C03Q coa/3 + Acy sina 5 tn/?} 
“—{—Ac* sina sin0 + Acy sina cos^} 
{— Ae* sina cos/3 + Acy coaa sin/Sy 


(2.143) 

(2.144) 

(2.145) 


The right hand size of the system of first order equations (Eqs. 2.36-2.38) consists 
of harmonic forcing functions ho, hr, and These functions are 

essentially harmonic functions from Eqs. (2.143-2.145) Based on this fact, the solution 
is assumed to be harmonic functions of a and /?. 


Plx 

= ai{z,/3)cosa 

+ 

a2{z,/3)sina 

(2.146) 

Ui, 

= a3{z,/3)cosa 

+ 

a4{z,/3)sina 

(2.147) 

Viz 

= ai{z,/3)cosa 

+ 

oo{z,P)sina 

(2.148) 

Ply 

= bi{z,^)cosa 

+ 

b2{z,P)sina 

(2.149) 

Uly 

= b3{z,/3)cosa 

+ 

b^{z,/3)sina 

(2.150) 


= bi[z,0)cosa 

+ 

bo{z,P)sina 

(2.151) 


Using the above substitutions in the set of first order equations Eqs. (2.36-2.38) 
yields 12 coupled linear partial differential equations. This set of equations are given 
in Appendix A, for the constant properties model. 

The solution procedure for the 12 linear PDE’s is exactly the same as that of 
the zeroth order solution. The integration is performed with a 4-5th order Runge- 
Kutta method, predictor-corrector method and Adams methods. All the methods 
yield almost identical results, with the Runge-Kutta based method being the fastest. 
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2.9.1 Comparison with Nelson and Nguyen’s Approach 

The original analysis assumed the variables, and bi to be harmonic and separated 
them into two auxiliary functions of the form, 

‘*•(^5^) = fi{z)co80 + gi{z)sin!3 (2.152) 

where fi[z) and gi[z) are assumed not to vary with Nelson and Nguyen (1988a, 
1988b) thereby apply a second separation of variables substitution to the first order 
differential equations. While the above form of assumed solution yields results that 
agree with available experimental results, an examination of the numerical values of 
the functions /,(z) and p,(z) revealed a /? dependence, particularly at eccentricities 
above 0.5. The inclusion of the circumferential gradients of these variable should 
therefore improve the solution at higher eccentricities. 

The a,- and 6, in the current analysis are totally general functions of z and /3 
which thereby avoids the mathematical contradiction discussed above. Furthermore, 
in many cases the results of the current approach show better agreement with exper- 
unental results than the earlier results. 


2.9.2 Boundary Conditions of Assumed Variables 

The first order boundary conditions are expressed in the assumed solution variables 
are, 


ai(0,/?) 

0,(0,^) 

‘» 6 ( 0 ,/ 3 ) 

M0,/9) 


“(1 + ^«) A >® 3 ( 0 ,^) 
-(1 + ^.)/> a 4 ( 0 ,/?) 


(2.153) 

(2.154) 

(2.155) 

(2.156) 


0 
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oi(£,/3) 

= 0 


(2.157) 


= 0 


(2.158) 

h{0,l3) 

= -(1 

-h6)p6s(0,/3) 

(2.159) 

6j(0,^) 

= -(1 


(2.160) 


= 0 


(2.161) 


= 0 


(2.162) 


= 0 


(2.163) 


= 0 


(2.164) 


2.9.3 Solution of First Order Equations 

The same solution procedure that is used for the zeroth order equations is used to 
solve the reduced first order equations given in Appendix A subject to the boundary 
conditions of Eqs. (2.153-2.164) for variables o,(z,/3) and bi{z,/3). 

2.10 Determination of Dyn aTni r Coefficients 

In this section, dynamic coefficients are derived form first order pressure distribution 
Pi(z,/3,t). The following linearized force-motion model for a 2-DOF vibration is 
used to define the rotordynamic coefficients. In this equation, Ax and Ay define the 
displacement of the rotor relative to a static operating point and AF*, AF„ are the 
components of the perturbed force due to first order pressure field, pi{z,P,t) The 
significance of each of the linearized coefficients have been explained in Chapter I. 


AF. 

AF„ 


> = 


— Ay* Kyy 



+ 


^XX Pxy 

— Cy* Cyy 
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r 1 


' \ 



i 

A£ 1 

-TTly* 

TTlyy 


Ay J 


(2.165) 


The perturbed or incremental force components acting on the rotor due to a 

small motion about a static eccentric position (i,y) is given by integrating the first 
order pressure field, 


2ir 

epicos0 Rd0dz (2.166) 

2ir 

episin/3 Rdfidz (2.167) 

The perturbation motion is described earlier by an elliptical orbit. The dis- 
placements, velocities and accelerations at any point on this elliptical orbit are given 

ty, 




Ai = Xcosu/t 

(2.168) 


Ay = Y simvt 

(2.169) 


Ax = —wXsintjjt 

(2.170) 


Ay = —wYcosu;t 

(2.171) 


Ax = —iJ^Xcosuit 

(2.172) 


Ay = —w^Ysivwt 

(2.173) 

At ujt = 0, sinuit = 

0, cosujt = 1 and Ay = Ai = Ay = 0. Substituting these 

values in Eq.(2.165), 

AF. 


A6s{ a** — } -f A6y{c^yO?} 

(2.174) 

AF. 


- A€,{-A:^ -I- -1- Ae„{C,^u>} 

(2.175) 


At wt - f, sinut = 1, cosuit = 0 and A® = Ay = Ai = 0. Substituting these values 



“t" 

Ae.{cy.u;} + Acy{K„ - My^u;^} 


(2.176) 

(2.177) 


AF, 

c. 

AFy 

c, 

R-om the following relationship between the first order pressure and the assumed 
variables, a,- and fc,-, 

epi — Ae,.pi, + Acypiy 

Pix and Ply may be expressed as, 

Pix = Oicosa -t- a2sina 

Ply = 6i C03Q + b2sina 

Equating Eqs.(2.166-2.167 and Eqs.(2.174-2.177) and dropping the perturba- 
tions, the following relations for the dynsunic coefficients are obtained. 


Kgx ^ ~ J J (^icos0Rd0dz 

(2.181) 

Cxyu; — — / / hicos/SRd/Sdz 

c» Jo Jo 

(2.182) 

5 1 

~ ^ "I" ~ ^ /o J 

(2.183) 

CyyU = — 1 bisin^Rdfidz 

Jo Jo 

(2.184) 

— Cxx<^ ~ ^ Jo Jo 

(2.185) 

L 2 1 

kyx — rrixyU} — — hiCosfiRd^dz 

(2.186) 

1 fL ,2ir 

Cyapu; — / / OL^sin^RdBdz 

Jo Jo 

(2.187) 

Kyy — MyyUj^ = — ^ ^ h2sin(5 Rd(5 dz 

Jo JO 

(2.188) 


(2.178) 

(2.179) 

(2.180) 
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These 8 equations must be evaluated for at least two whirl frequencies to obtain 
solutions for the 12 dynamic coefficients. A least squares approach is employed for this 
step. Typically, 2-4 whirl frequencies are used in a least squares scheme to compute 
these coefficients. Also, the dynamic coefficients for an anmil uT geal are essentially 
independent of whirl frequencies. The 2-D integration performed numerically is an 
improvement over the average value approach employed by the Nelson sind Nguyen 
(1988a, 1988b). 

2.11 Dynamic Coefficients baised on Extemsd Load Specification 

TypicaDy, for seals the dynamic coefficients are computed as a function of eccentricity. 
This assumes that the eccentric position of the shaft has been specified and the 
resultant reactive force due to the pressure distribution in the seal is to be determined. 
In some cases, it is possible to specify the angle at which external load is supported 
by the seal during the operation of the turbomachine. For example, unit 3-01, an 
experimental seal under design at NASA (results to be discussed later) supports the 

external load at a constant angle of 290° in the rotor coordinate system as shown in 
Figure 2.15. 

The problem now is to determine the eccentric position given an external load 
F and its load angle #. This is accomplished by iteratively searching for an eccentric 
position e of the rotor which produces a pressure distribution po{z,/3) which when 
mtegrated over the entire seal balances the applied load in magnitude and direction. 
This iterative search is carried out using a modified 2-D Newton-Raphson method 
discussed below. 
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Figure 2.15 Ex a m ple of External Load acting on a Seal 
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Figure 2.16 External Load and External Load Angle 
2.11.1 Steady State Force Equilibrium Position 

A modified 2-D Newton-Rapbson method is used to locate the operating position. At 
the steady state equilibrium position, 

Let Fg and Fy be the X and Y components of the seal force obtained by inte- 
grating the pressure field within the seal for a given rotor position (x,y). Let F* and 
Fy be the components of the external load. 

/. = F. -I- F. (2.189) 

fy = Fy + Fy (2.190) 

where /, and fy are the residual forces in X and Y directions respectively. The 
problem reduces to finding (x,y) such that the residual forces /. and /„ are zero. In 
other words, find a rotor eccentric position such that F* is balanced by F, and Fy is 



balanced Fy . 

Tbe iterative search is described by, 
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^k+\ = ifc + Axk 

Pk+i = y* -f Ay* 


(2.191) 

(2.192) 


and the increments Ai* and Ay* are computed from Eq. (2.193). 


efk 

3x 

BFy 

6x 


By 
BFy 
^ . 


Az* 

Ay* 


Sv 


(2.193) 


The derivatives in Eq. (2.193) are computed using a forward difference formula as 


^ ^ i".r(x + Az,y) - F,(z,y) 
dx Az 

^ ^ i",(g,y + Ay) - F,(a,v) 
^y Ay 

^ ^ i=;(x + Az,y) - F„(z,y) 
dx Az 

^ ^ n(g.y + Ay) - FJz,y) 
dy Ay 


(2.194) 

(2.195) 

(2.196) 

(2.197) 


The iterative search stops when the residuad forces /* and fy lire below some 
specified tolerance. Once, the eccentric position is determined, the computation of 
d3mamic coefficients is carried out as before. 
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CHAPTER III 

VARIABLE PROPERTIES MODEL 

The effect of variable fluid properties as related to liquid seals for cryogenic applica- 
tions was first investigated in some detail by Simon and Rene (1989). Their initial 
work did not include fluid inertia effects and the analysis was based on a simplified 
Reynolds equation. San Andres (1991) developed a seal analysis that included vari- 
able fluid properties as a function of local pressure and a mean temperature. He 
used the NIST 12 Database, MIPROPS (1986) to compute the fluid properties. This 
database is based on the 32-term Modified Benedict- Webb-Ruben Equation of State. 
Experimental and theoretical data is used to compute the coeflRcients of the terms in 
this equation of state and data is available for a number of fluids. 

The working fluid in SSME turbopump is either liquid oxygen (LOX) or liquid 
hydrogen (LH2). Figure 3.1 shows the variation of density and viscosity of LOX a 
function of pressure. Typical inlet pressures for the turbopump are in the range of 
20 Mpa and the exit pressures are in the 3 Mpa range. 

® typical seal, inlet and exit conditions are given below. 


inlet pressure, p, 19.0 Mpa 

exit pressure, pt 3.0 Mpa 

mean temperature, T* 90° K 

For the above conditions, the fluid properties for LOX at inlet and exit are 


at inlet: 
density, p,- 
viscosity, /i,- 


1179 kg/m® K 
2.32x10-“ Pa-s 
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compressibility, 0.0014 1/MPa 

at exit-. 

density, p, 1148 kg/m^ K 

viscosity, p, 2.01x10"“ Pa-s 

compressibility, 0.0019 1/MPa 

For the same conditions, the fluid properties for LH2 at inlet and exit are, 
at inlet: 

density, pi 42.13 kg/m® K 

viscosity, p, 6.43x10"® Pa-s 

compressibility, 0.03335 1/MPa 

at exit: 

density, p. 8.18 kg/m® K 

viscosity, pe 4.13x10"® Pa-s 

compressibility, 0.334 1/MPa 

The change in density and viscosity for LOX is relatively small (about 2.6% 
and 12% respectively for the above case). However, for LH2, the change in density 
between inlet and exit is considerable and this will have a noticeable effect on the 

dynamic coefficients computed (San Andres, 1991). A similar change is also noticed 
in viscosity. 
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Figure 3.1 Properies of Liquid Oxygen at 90°K 
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PrtMur* (MPa) 



Figure 3.2 Properies of Liquid Hydrogen at 90°K 
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3.1 Thcnnopliysical Properties Model 

Even though the following analysis is valid for any liquid, the two fluids of interest in 
this research are the cryogenic fluids, liquid oxygen and liqxiid hydrogen. These are 
the two working fluids commonly used in the SSME turbopump. 

The standard model generally used for representing the thermophysical proper- 
ties of fluids is the "Modified Benedict-Wcbb-Ruben (MBWR)” equation of state. 
This particular model is widely used to correlate thermodynamic property data and 
a number of computer codes are available to tabulate the properties for various fluids 
based on this equation of state. The most important of these codes is NIST Standard 
Rtftvtnct Database 12^ published by National Institute of Standards and Technology. 
This code is available in source form and can be easily integrated into a seal code and 
this code in its source form is used in the present work. 

The main advantages of a MBWR based fluid property model are, 

1. Accurate data available for a number of fluids. 

2. Easy adaptability to use in a seal code. 

3. Correlation of experimental data from various sources. 

4. Lends itself to analytical work. 

3.2 MBWR Equation of State 

The MBWR equation of state in the single phase region is a 3 2- term equation given 
below. 


P 


pRT + 

P^{G{l)T + G(2)r^/* + G{3) + G(4)/T + G{5)/T^) + 
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^®(G(6)r + G{7) + G{8)/T + G{9)fT^) + 

/(C?(10)r + G(ll) 4- G{12)IT) + /(C?(13)) + 

/(G(14)/r + G(15)/r*) + /(G(16)/T) + 

/)*(C?(17)/r + G(18)r*) + p\G{ 19 )/T^) + 

P^{G{20)IT^ + G(21)/r®) +. 

p^{G{ 22 )/T^ + G(23)/r^) + 

/(G(24)/r^ + G(25)/r») + 

p®(G(26)/r* + G{27)IT*) + 

p“(G{28)/r* + G(29)/T®) + 

p'^(G(30)/T* +. G(31)/T'* + G{Z2)IT*) (3.1) 

where, 

p pressure 

p density 

T absolute temperature 

7 Pc density at Tcrit 

^(*)>* = li2, ...32 linear coefficients 


The li n ea r coefficients G^(t) are computed using experimental and analytical data 
for various fluids. 

The expression for viscosity is given as, 

= Mo(T) + pi{T)p + P2{p,T) 

Mo = S G4i)T("-')/* 

t=l 

m{T) = /’.(l) + F.(2) {F.(3) - /n(T/F.(4))}’ 


(3.2) 

(3.3) 

(3.4) 
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F(p,T) = £,(1) + E^(2)H(p) + E,(3V"‘ + 

£.(4)H(p)/r= + B„(5)/Vr‘‘ + 
B„(6)/T + E,{7)H(i>)IT 
G(T) = £^1) + B„(2)/T 
B(p) = p'‘‘(p - £48))/£„(8) 


(3.5) 


(3.6) 

(3.7) 

(3.8) 


3.3 Bulk Flow Governing Equations 


The bulk flow governing equations for compressible flow are given by Nelson (1985). 
The original equations are derived for gas seals and the same equations will be used 
for this analysis using MBWR equation of state. 

Continuity: 


d{phu) 1 d{phv) d{ph) _ 

dz Rd/3 dt ~ 


(3.9) 


Axial Momentum: 


h dp 
p dz 


,fdu V du 
+ ft-zVu^ + + 



(3.10) 


Circumferential Momentum: 


pRdp ~ 

+ +v2 + /r ~ + (t; - w)^ (3.11) 


A coxnpanson of the above goveming equations with the case of constant prop- 
erties model (£qs. 2.13—2,15) reveals that these goveming equations are essentially 
the same except for the continuity equation where the density term is retained within 
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the derivative. The momentum equations are the same, but with a variable density 
and viscosity. 

3.4 Comparison with San Andres (1991) 

San Andres (1991), presented an analysis for variable properties based on a finite 
difference formulation. The solution is based on a finite difference scheme that is based 
on a method of Launder and Leschziner (1978) and used a SIMPLEC algorithm of 
Van Doormal and Raithby (1984). In the current work, a completely different solution 
procedure will be used. The analysis developed for the constant properties case in 
the previous chapter will be extended to the case of variable properties. 

3.5 Perturbation Analysis 

In this section, the set of governing equations, Eqs. (3.9-3.11) are perturbed about 
their steady state values to obtain the zeroth and first order equations. The procedure 
is similar to the one outlined in Section 2.5. 

The assumed form for the dependent variables, film thickness and the fluid prop- 


erties for perturbation are given as, 

= uo(z,/3) + €Ui(z,/3,0 (3.12) 

v{z,!3,t) = wo(^,/3) + evi(z,/9,t) (3.13) 

p(«,/3,<) = Po(z,/?) + epi(z,/3,<) (3.14) 

^(^»^»^) = + ehi(z,0,t) (3.15) 

p{z,/3,t) = po{z,l3) + cpi(z,/3,t) (3.16) 

= fio{z,/3) + (3.17) 


where tto, vo» Po, ho, po, Po are the zeroth order variables and «i, vi, pi, hi, pi, pi 
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are the corresponding first order variables, and eui, eri, epi, efij, epi, tpi are the per- 
turbations. The zeroth order variables, po tmd po and the corresponding first order 
variables pi, pi are not independent and will be related to the primary variables u, 
V, p through the MBWR equation of state. Substitution of these expressions into 
Eqs. (3.9-3.11) and neglecting second and higher order terms }rields the sets of zeroth 
order and first order equations. 

(zeroth order equations) -f e(first order equations) = 0 


3.5.1 Zeroth Order Equations 

The zeroth order equations are given by. 
Continuity. 


d{pohoVo ) 1 d{pohoVo) /a io\ 

dz R dp ~ 

Axial Momentum: 


Hq dpo 
Po dz 


ho{ 


Vo duo 


+ 


5uo, 


R dp 

+ + ^0 + fro-^y/ul (vo — tv)* 


(3.19) 


Circumferential Momentum: 


ho dpo 
PoR dp 


Vo dvo 




dvo, 
« 0 — } 


^°^Rdp ' "“dz 

+ /»0Y\/“0 + + frO ^ ° \/t^o (^0 - tv)^ 


(3.20) 


2 
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3.5.2 First Order Equations 

The first order equations are given below. 
Continuity: 


PqVq dhi poho dvi voho dp 


R d(3 


+ 


R dfi 


+ 


R dp 


dui 

dz 


dp 

+ />oho-^— + lioho-^ 


dz 

dh 

-\-A^ui + ApPi = — />o“o~5 — “ Po~^ ~ Ahhi 


dhi 


dz 


dt 


(3.21) 


Axial Momentum: 


hovo dui ^ dui , dpi , dui „ _ _ 

+ = Bfchi (3.22) 


Circumferential Momentum: 


JiqVo dvi ho dp 


+ 


R dp po dp 


, dvi dvi 
+ houo-^ho-^ 


+ CuUi + C'vUj + C/,p\ 

+ Cf^Pi = Cfihi (3.23) 


where the coef&cients A^, Ap, • • • etc. are functions of steady state variables uo, Uq, 
Po and their axial £ind circumferential gradients. 


3.6 Zeroth Order Boundary Conditions 

The boundary conditions are similar to the constant properties model, except for 
exit pressure term which is retained for the variable properties case, since the fluid 
properties vary with pressure. The boundary conditions for zeroth order equations 
are shown in Figure 2.6. 
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The following is a summary of the boimdary conditions for the zeroth order 
equations. 

At the inlet: 
axial velocity, uq: 
prior to inlet: 

uo(0,/?) = 0 ' (3.24) 


right after inlet: 


^o(0,/?) — Uoi(0,/3) 


(3.25) 


circumferential velocity, uq: 


uo(0,/3) = psrxuiR (3.26) 

pressure, p^: 
prior to inlet: 

Po(0,/3) = Pi (3.27) 

right after inlet: 

Po(0,;3) = poi(0,^) (3.28) 

The pressures pi, Poi(0,/?) and axial velocity uoi(0j/?) the inlet are related by, 


1 


Pi - Poi(0,/3) = rM0,/3)«2i(0,^)(l + e0 


(3.29) 


At the exit, the exit pressure recovery coefficient is assumed to be 1, i.e., 


Po 2 ( 0 ,^) - p, = rpo(0,/3)u^,(0,;9)(l-l) 


(3.30) 


or, 


Po2(0,/3) = pt 


(3.31) 
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Eq. (3.29) may be revrritten as, 

Poi(0,/3) = Pi - ^/>o(0,)3X(0,^)(l + ^i) (3.32) 

or 

>^.(0,^) = + (3-33) 

At the outset, poi(0,/3) is unknown and must be solved iteratively by requiring 
that the pressure distribution at the seal exit satisfies the following condition. 

Po2{L,/3) = pe (3.34) 

subject to the constraints of Eqs. (3.26,3.29) 

3.7 Reduction of Zeroth Order Equations 

In the following analysis, the original zeroth order equations are reduced into a form 
suitable for the solution procedure developed in Chapter II. 

The fluid properties are, in general, functions of local pressure and temperature 
as given below. 


P = p(p,T) (3.35) 

p = p(p,T) (3.36) 

The dependent variables, pressure p and temperature T are functions of the a-rUI 
and circumferential coordinates, («,/3). 


P = 


T 


P(^,/?) 


(3.37) 

(3.38) 
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or Eqs. (3.35,3.36) may be rewritten as, 

P = pipiz,fl),T{z,P)) (3.39) 

P = fi{p{z,0),T{z,/3)) (3.40) 

Using chain rule for differentiation, the terms which represent 

the changes in fluid properties with respect to axial and circumferential coordinates, 
are expressed (subscript 0 refers to zeroth order variables) as, 


dpo 

_ dpo dpo dpo dTo 


dz 

dpo dz * dTo dz 

(3.41) 

dpo 

_ dpo dpo dpo dTo 


d0 

dpo d(3 ^ dTo d/3 

(3.42) 

dpo 

_ dpo dpo dfio dTo 


dz 

dpo dz dTo dz 

(3.43) 

dpa 

_ dpo dpo dpo dTo 


dp 

dpo dfi dTo dfi 

(3.44) 

In the following einalysis. 

fluid properties are assumed to be 

a function of the 


local pressure and a mean temperature, T*. 


P = P(p,r*) (3.45) 

P = M(p,r*) (3.46) 


In other words, the fluid flow is treated as an isothermal flow and the flmd 
properties, density and viscosity are are assumed to vary as function of local pressure, 
p(^»/^)j oidy. For a constant temperature field, the partial derivatives with respect to 
temperature T in Eqs. (3.41—3.44) vanish giving the following simplified expressions. 


dpo 

_ dpo dpo 

dz 

dpo dz 

dpo 

_ dpo dpo 


dpo dfi 


(3.47) 

(3.48) 
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_ dpo dpo 

dz 

dpo dz 

dpo 

_ dpo dpQ 

d/3 

dpo dfi 


(3.49) 

(3.50) 


Sfit 

®Po 

§JiS. 

Bpo 

ax 

§21 

aa 


RaIc of cHaixge of density with pressure 
R&te of change of viscosity with pressure 
Axial pressure gradient 
Circumferential pressure gradient 


The term ^ is related to the compressibility of the fluid and is usually represented by 
the dimensionless parameter, isothermal compressibility, ^^|r- A larger isothermal 
compressibility signifies a more compressible fluid. 


Using these relations, the Eqs. (3.18-3.20) may be rewritten as given in Appendix 
B, Eqs. (B.1-B.4). For constant properties, i.e., an incompressible fluid, po,Mo = 
constant, or, ^ ^ = 0, and the Eqs. (B.1-B.4) reduce to the Eqs. (A.1-A.4) of 

the constant properties model. 


As in the case of constant properties model, the reduced zeroth order equations 
may be rewritten with all /^-dependent terms on the right-hand side as, 

K(uo,vo,po, ^) 

(3.51) 

Fp(uo, Wo.Po, 1^) 

The fimctions F„, Fp for the variable properties model are given in Appendix B. 



3.8 Solution Procedure for Zeroth Order Equations 

The solution procedure for zeroth order solution is exactly the same as discussed in 
Section 2.7. The only difference is that the fluid properties and their dependent terms 
are updated at each grid point during numerical integration. 
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3.9 Reduction of First Order Equations 

The original first order equations given in Eqs. (3.21-3.23) are reduced the 

zeroth order equations. 

The first order variables in the first order equations, Eqs.(3.21-3.23) are ui, vi, 
P\i Pit Out of these variables only Ui, vi, pi are primary variables. The remaining 
two variables pi, are related to pi using the property relations based on the MBWR 
equation of state. 

The relationship first order variables pi and p\ is given by 


Pi — ^piPi 


(3.52) 


where, 


^ 


(3.53) 


The above relationship is obtained by perturbing the MBWR equation of state with 
respect to pressure and density and equating the terms on both sides. Prom a different 
perspective the ratio ^ is a ratio of two infinitesimally small quantities, which is 
nothing but the derivative 1^, 

&p 

In the expression for viscosity given in Eqs. (3.2-3.S), the relation between pres- 
sure and viscosity is not explicit. Let the first order variables pi and pi be related 


where. 


Pi — Sftipi 


R _ 

Rui = ^ — 


(3.54) 


(3.55) 


The first order variables pi and are then related through. 


Pi = 


(3.56) 
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where, 

(3.57) 

The terms ^ £uid are obtained by differentiating Eqs. (3.1,3.2-3.8) with re- 
spect to density p and viscosity p respectively and these expressions are given in 
Appendix F. 

Using the relations Eqs. (3.52-3.57), the first order equations are rewritten in 
terms of the primary first order variables ui, vi and pi. 

Continuity: 


PqVq dhi poho dvi , voho dpi 


R d(3 


+ 


R d!3 


R d0 


+ Poho 
dhi 


dui . , 

87 + + 

dhi 


hoApi 


+A^ui + A^vi + ApPi = -pouo^ po-^ Ahhi 


dz 


dt 


dpi 

dt 

(3.58) 


Axial Momentum: 


hoVo dui 

R dp 


. dui 


hp dpi 
Po dz 


+ ho 


dui 

~dt 


-h -t- BvVi -F Bj,pi 


= Bnhi (3.59) 


Circumferential Momentum: 


hpVo dvi 

~Wd0 


hp dpi 

ppR d0 




dvi 


CftUi Cx)Vi -|- CpPi 


= Chhi (3.60) 


The coefficients A„, Ap.. for a variable properties model are defined in Appendix B. 
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3.10 First Order Boundary Conditions 

Perturbing tbe zeroth order boundary conditions of Eqs. (3.26,3.29,3.34) yield the 
following first order boundary conditions. 

At the inlet: 

ctrcumfereniial velocity, vi: 


vi{0,(3) = 0 


pressure, p^: 


(3.61) 


j{2/>o(0,/3)a„(0,/3)u,(0,/3) + u5,(0,^)p,(0,/3)) (3.62) 

At tbe exit: 
pressure, pi : 


Pi(0,^) = 0 (353) 

Using the relation between pj and pi, Eq. (3.62) may be rewritten as, 

{1 + 0.5(1 + ^,)ugi(0,)3)^}^^^ ^‘)P“(®’^)“o(0i/?)}«i(0,/3) (3.64) 

3.11 Solution of First Order Equations 

The same procedure developed for the constant properties model is used for this case. 
The set of first order equations in Eqs.(3.58-3.60) are further reduced by employing 

“separation of variables” technique for an assumed small motion of the vibrating 
rotor. 

The right hand rize of the system of first order eqnations consists of harmonic 
forcing fnnctions f.. Based on this fact, the sointion is assumed 
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to be harmonic functions of a and /3. 


Plx 

= ai(z,0)cosa 

+ 

a2{z,0)sina 

(3.65) 


= a3(z,/3)cosa 

+ 

a4{z,(3)sina 

(3.66) 


= as(z,/3)cosa 


a«{z^0)3ina 

(3.67) 

Ply 

= bi(z,/3)cosa 

+ 

b2{zy0)sina 

(3.68) 

«lv 

= b3(z,j3)cosa 

+ 

b^(z,fi)aina 

(3.69) 

Vly 

= b3[z,f3)coaa 

+ 

6«(z,/?)stna 

(3.70) 


Using the above substitutions in the set of first order equations Eqs. (3.58-3.60) 
yields 12 coupled linear partial differential equations. 

The first order boundary conditions expressed in the assumed solution variables are: 


ai(0,/3) = 


a,(0,/?) 

«b(0,/3) 

oe(0,/9) 


= 0 
= 0 
= 0 


(l + 6)po(0,yg)a3(0,^) 

{1 + 0.5(1 + e0<(0,^)g} 

{1 + 0.5(1 + e,)ug,(0,^)g} 


a,(L,/3) 

ii(0,/3) 

62(0, / 9 ) 
65(0,^) 
be(0,/3) 


0 

{1 + 0.5(1 + 6)ug,(0,/?)g} 
{1 + 0.5(1 + ^4«2,(0,y9)g} 

0 

0 


(3.71) 

(3.72) 

(3.73) 

(3.74) 

(3.75) 

(3.76) 

(3.77) 

(3.78) 

(3.79) 

(3.80) 
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6i(I,/3) 
62(1, / 3 ) 


( 3 . 81 ) 

( 3 . 82 ) 


The solution procedure follows the same steps as explained in Chapter II, 
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CHAPTER IV 

THERMAL EFFECTS MODEL 

Yang ti al. (1992) and San Andres et al. (1992) conducted a thorough investigation 
of thermohydrodynamic (THD) analysis for cryogenic seals. Yang developed an ap- 
proximate THD analysis and provided steady state solution for the case of a centered 
seal. San Andres (1992) introduced a full set of bulk flow governing equations for 
THD imalysis and investigated it with a finite difference based numerical solution. 
The results from their study show that for some cryogenic fluids temperature rise due 
to friction at the rotor surface may lead to a two phase flow, seriously affecting the 
performance of the turbomachine. 

The goal of this work is to extend the solution procedure developed in previous 
chapters to a THD analysis. The governing equations used for this suialysis are based 
on San Andres et al. (1992) and are given in Eqs. (4.1-4.4). In the analysis developed 
in this chapter, the zeroth order equations will be solved for the primary variables, 
Uoy Vo, po and To for a centered seal. The temperature distribution will then be used 
with the variable properties model developed in Chapter III to compute the d}mamic 
coeflicients. In other words, pertiirbation due to temperature are not included in the 
analysis. The goal is to show the viability of current analysis to handle thermal effects 
of the energy equation. 

Continuity: 

Axial Momentum: 

d , 9 , 2\ ^ / I \ ih 

^(puk) + j-^(phu ) + = -kgj + r«|o 


(4.2) 
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Circwmfertntial Momentum: 


d , , \ 5 / , X ^ / , 5\ i^P ih 

-^(Khv) + ^^(phuv) + + T„lo 


(4.3) 


Energy Equation: 


CA§i{phT)^.l-^pkuT) + .l-^(phvT)} + g, = + 

+ - v(r„|J) - ti(r.,,|J) (4.4) 


where, q = R/3 and the shear stress terms, r^, etc., are given as. 


^*vlo = -/>(/. + V* + /r|\/«^ + (w-tw)2) 

= + /r — ^ yju^ + {v- u;)^) 


'9VlO 


r 1. - _ / (^ - ^) /? 1 

“ 2 a? ^ 4h^’''2^' 2 


The heat transferred, Q, is given by. 


C?. = 0 

= hr{T — Trotor) + h,{T — T^al) 


(4.5) 

(4.6) 

(4.7) 


(4.8) 

(4.9) 


Eq. (4.8) is for adiabatic case where there in no heat transfer, while Eq. (4.9) is for a 
general case. The stator and rotor Reynolds numbers are given as. 


R, = + v’ (4.10) 

Rr = + (tJ — U>y (4-11) 

The heat transfer coefEcients for rotor and stator are given by, 

h.{z,q) = (4.12) 

KM = pC,^u‘ + (v-wr[^)(f!^)-i (4.13) 
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The fluid properties, density p find viscosity p are assumed to be a function of 
local pressure and temperature. MWBR equation of state is used to represent the 
fluid properties. 

The variation in fluid properties with respect to the eudal and circumferential 
coordinates, ^ are expressed as (Eq. 3.39-3.40), 


dpo 

dz 

dpo 

dl3 

dpo 

dz 

dpo 

dH 


dpo dpo 
dpo dz 
dpo dpo 

dpo d(3 

dpodpo 
dpo dz 
dpodpo 
dpo d(3 


dpo 

dTo 

dTo 

dz 

dpo 

dTo 

dTo 

dP 

dpo 

dTo 

dTo 

dz 

dpo 

dTo 

dTo 

dfi 


(4.14) 

(4.15) 

(4.16) 

(4.17) 


where *ire the circumferential gradients of the primary variables, 

u,v,p, and T respectively. 

Using these relations, Eqs. (4. 1-4.4) may be rewritten as. 


dh , dh dh ^du \dv dT,^dp^ dT,^ dp. 


"'‘S + "'‘"SI = "'■S + 

'’'■I + = -'‘I + 


„.,8T dT BT^ . 


+ fa;R(r^|h) — v(t^Io) ~ '*(‘Trvlo) (^-^l) 


where, 

p density 
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i£ 

Bp 

§£. 

BT 


Pv 

k 

f. 


fr 


dynamic viscosity 

rate of change of density with pressure 

rate of change of density with absolute temperature 

specific heat at constant pressure, |^|p 

volumetric expansion coefficient, 

thermal conductivity of fluid 

stator friction factor 

rotor friction factor 


4.1 Zeroth Order Equations 

The zeroth order equations are obtained by dropping the time dependent terms in 
Eqs. (4.18-4.21). 

The above equations may again be written in the following fashion. 


dz 


Fu(tto, vo,po, To, 

&vo 

dx 


Fr(Uo, Uo,POt To, ) 

^Po 

Bz 


T'p{uo,vo,po,To, 

BTti 

Bz 


Fr(tio,wo,Po,ro, 1^, 


The functions F„, Fp, Ft are given in Appendix C. 

4.2 2^oth Order Boundary Conditions 

The boundary conditions for the zeroth order or steady state equations is illustrated 
in Figure 4.1. 

The following is a summary of the boundary conditions for the zeroth order 


equations. 
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Figtire 4.1 Zeroth Order Boundary Conditions 


At the inlet: 
axial velocity, Uo: 
prior to inlet: 

right after inlet: 

circumferential velocity, vq: 


uo(0,/?) = 0 


uo(0,/?) = uoi(0,/3) 


(4.23) 


(4.24) 


Uo(0,/?) = psrxu/R 


(4.25) 


pressure, p^: 
prior to inlet: 


Po{0,p) = Pi 


(4.26) 
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right after inlet: 


Po(0,/3) = poi(0,/3) 

temperature, Tq: 
prior to inlet: 


(4.27) 


3o(0,^) = T, ( 4 . 28 ) 

right after inlet: 

Jo(0,/?) = Toi (4.29) 

The pressures pi, poi{0,(3) and axial velocity uoi(0,/?) at the inlet are related by, 


Pi - Poi(0,/3) = 2Po(0,)9)u2i(0,/?)(l + ei) 


(4.30) 


To(0,/?) 


Ti- 


2Cp. 


(1 + Ci){l - 


Po{0,/3) 

Pi 


(1 - Ti/3„,)} 


4.3 Solution Procedure for Zeroth Order Equations 


(4.31) 


The solution is based on an iterative procedure and the steps involved are explained 
below. 

The continuity and momentum equations are solved using the solution procedure 
developed in previous chapters. At the outset, a nominal temperature distribution is 
assumed for the entire flow field. TypicaUy, a constant temperature, temperature 
at inlet, is assumed. At the end of convergence of the continuity and momentum 
equations for a given temperature field, the energy equation is integrated using the 
updated variables tt, v and p. from the previous iteration. At the end of one complete 
integration of energy equation, a new temperature distribution, T^{z,(3), is available 
and it is used to update the fluid properties. 

The set of contmuity and momentum equations are again solved with the up- 
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dated properties until the solution converges. The new tt, r, and p are then used in 
the integration of the energy equation to obtain a new updated temperature distribu- 
tion. The cyclic procedure is repeated until two successive temperature distributions 
converge to a specified tolersince. 

Typically, it takes about 6—8 iterations for the solution to converge. 

4.4 Comparison of Current Analysis with San Andres et al. 

San Andres et al. (1992), uses a finite difference scheme to solve the coupled, nonlinear 
PDEs of continuity, momentum and energy equations. The procedure they use is 
based on the forward marching solution of Launder and Leschziner (1978) and uses 
SIMPLEC algorithm of Van Doormal and Raithby (1984). The flow domain is divided 
into control volumes and governing equations are integrated on the control volumes to 
give sets of no nlinea r algebraic difference equations for each primary vzu'iable. Then 
they implement an iterative procedure where the continuity and momentum equations 
are solved followed by the energy equation. The continuity and momentum equations 
are allowed to converge to an intermediate limit and then the energy equation is 

solved. Fluid properties are updated and the procedure is repeated until the solution 
converges. 

The main advantage of the solution procedure used in the current work is its 
simplicity. 

In the current work, an iterative procedure based on the 
direct integration of governing equations, is implemented. The problem is reduced to 
repeatedly solving a system of 3 ODEs until the temperature distribution converges. 

San Andres et al. reported that their method requires about 20 iterations for 
the solution to converge. Current solution procedure takes about 6-8 iterations to 
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converge. 

4.5 First Order Solution 

The equations used for first order solution are the same as variables properties model. 
As mentioned earlier, the governing equations are not perturbed with respect to 
temperature T . The temperature distribution obt^ed from the zeroth order solution 
is used to update the properties as a function of local pressure and temperature and 

the dynamic coefBcients are computed using the variable properties model of Chapter 

III. 
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CHAPTER V 

ARBITRARY PROFILE SEALS 

The typical seal geometry of an interstage seal of SSME turbopump has either a 
straight or a tapered (convergent) axial profile as shown in Figure. 1.2. This nominal 
seal profile may be altered during the course of its operation, for example, as in 
the case of SSME interstage seals due to mechanical and thermal distortions. An 
example of predicted seal profile of an interstage seal is shown in Figure 1.3. Tests 
at NASA/MSFC reveal that seals initially designed with a large convergent taper 
have become divergent over a part of the length of the seal during the course of their 
operation changing the dynamic characteristics of that seals. 

For the purpose of this study, an arbitrary profile seal is a seal whose geometry in 
axial and circumferential directions may vary in any specified fashion. For example, 
the distorted seal shown in Figure 1.3 would be an example of an arbitrary profile 
seal. 

This change in the seal profile may be due to distortion such as in the interstage 
seals of SSME turbopump or by design itself in order to enhance some optimum 
dynamic characteristics of the seal. It has been known for a long time that convergent 



zeroth degree first degree second degree third degree 


Figure 5.1 Examples of Seal Axial Profiles 
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seals provide a higher stiffness than a straight seal. It is ako a fact that divergent 
seals under certain conditions lose stiffness with eccentricity. Generally, of sdl the 
variables that go into the design of zin annular seal, the variable that can be easily 
modified is the profile of the seal whether it be a straight, tapered or more exotic 
shapes. As will be shown with reference to an elliptical seal, even a small change in 
seal profile can have a noticeable effect on the d}mamic characteristics of the seal. 
The fact that dynsunic coefficients can be varied with the profile may be used as a 
design criterion to build a seal with a set of optimum dynamic characteristics. 


5.1 £x2imple of an Arbitrs^y Profile Seal: Elliptical Seal 


Annular seals initially designed with a circular cross section have been found to ac- 
quire, under load, an oval shape similar to an ellipse. Currently, investigations are 
being carried out to study the effect of such a change in profile on the flow rates and 
dynamic characteristics of these seals. This is the motivation behind the following 
study of an elliptical seal. 

An elhptical seal is am annular seal with an elliptical cross-section as shown in 
Figure 5.3. The clearance function for this seal varies in the circumferential direction 
as opposed to a constant clearance for a typical straight or tapered seals with circular 
cross-sections. The degree of curvature of the elliptical seal compared to a circular 
seal is specified by the psirameter ellipticity, 6 as, 


e, -Cy 

c. 


(5.1) 


where c* and Cy are radial clearances at sexxu-niajor and semi-minor axes respectively. 

For ^ = 0, the seal is a circular seal and for 5=1 the stator contacts the rotor 
and for any value in between, i.e., 0 < 5 < 1, the seal is an elliptical seal. Figure 5.4 
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illustrates an elliptical seiil with various ellipticity factors. 

The following different cases of axial profile are studied in detail for this elliptical 

seal. 

1. Straight Profile 

2. Linear Profile 

3. Quadratic Profile 

These axial profiles correspond to zeroth, first order and second curves shown in 
Figure 5.1. 

The equation of sin ellipse is given by, 

X = acosfi (5.2) 

y = hsinfi (5.3) 

where a and h are the semi-major and semi-minor axes respectively. At any ^tTi giilnr 
position along the circumference, the radius r of the ellipse is given by, 

»•(«>/?) = \j{acosfif -I- {hainPf (5.4) 

and the clearance c at this location is given by, 

c(r,/9) = r(z,/3) - R (5.5) 

where R is the radius of the rotor. If the semi-major and semi-minor axes of the ellipse 
vary in some functional form along the length of the seal, the clearance function of 
this seal is given by. 




(5.6) 
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where fi{z) and fziz) are the semi-major and semi-minor axes variations along the 
z-axis of the ellipse. The function fi{z) for straight, linear and quadratic axial profiles 
are given below. 


/l(^) = 

fll 

(5.7) 

/l(^) = 

Oi -f U2Z 

(5.8) 

Mz) = 

Oi -f 022 + 032* 

(5.9) 

f2{z) = 

bi 

(5.10) 

f2{z) = 

bi -i- 62^ 

(5.11) 

f2{z) = 

^1 + ^22 -|- bsZ* 

(5.12) 


where Oj, oj, 03, 61, 621 are constants that define the curvature of the axial profile. 

The curvature of the elhpse in the circumferential direction is varied using the 
parameter, elhpticity S. The elhpticity of the seal is defined as (Figure 5.3), c* = c,- 
at inlet and = c, at exit. Therefore, the clearance Cy is given by, 

cy = c.(l-^) (5.13) 

The Appendix E provides the functions fi{z) and /2(x) for straight, liTieAr and 
quadratic axial profiles as a function of S. 

5.2 Results 

Three different axial profiles are considered for the elliptical seal and their dynamic 
characteristics are studied with reference to similar seals with circular cross-sections. 
Results are provided for two cases. 


1. Straight elliptical seal compsued to a straight circular seal. 
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2. Elliptical seal with a linear iixial profile compared to a similar seal with a curved 
(quadratic) profile. 

5.2.1 Straight Elliptical Seal vs. Straight Circular Seal 

The plots shown in Figures 5. 4-5. 7 are comparisons for a straight seal with two 
elliptidty factors. For, 5 = 0, the seal is a circular (straight). For, S = 0.4, the seal 
is an elliptical (straight) seal. The various coefficients are normalized with respect to 
straight circular seal (5 = 0). 

Figure 5.4 shows the variation of direct stiffiiess for both seals as a function of 
eccentricity ratio. Even though, both seals have roughly similar stiffness to start 
with, for straight elhptical seal, both and Kyy decrease with eccentricity. In 
other words, there is a loss in stif&iess with eccentric operation. 

The cross coupled sti&ess in Figure 5.5 shows the opposite trend, i.e, they 
mcrease with eccentricity, almost exponentially. Such a large in cross coupled coeffi- 
cients should be a cause for concern form a stability point of view. 

Damping, shown in Figure 5.6, increases slightly with eccentricity but not to the 
extent of the circular seal. Leakage, shown in Figure 5.7, is smaller for elliptical seal. 

In short, the elliptical seal compared with the circular seal shows the following 
trends as a function of eccentricity. 

1. Loss in direct stiffness. 

2. Large increase in cross coupled stiffness. 

3. Relatively small increase in damping. 

4. Reduction in leakage. 

Except for the reduced leakage rate, all other comparisons point to the fact that 




Eccentricity Ratio 


Figure 5.5 Normalized Cross Coupled Stiffness for Elliptical Seal (Straight), (^=0 
5=0.4) 
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a straight elliptical seal is a bad design compared to a similRr circular seal &om a 
rotordynamic point of view. 

5.2.2 Linear Profile vs. Curved Profile 

The results shown in Figures 5.8-5.11, compare the results of an elliptical seal with 
a linear axial profile with a similar with a curved (quadratic) axial profile. The 
results are presented for a centered seal as a function of elliptidty, S. The dynzunic 
coefficients are normalized with respect to the coefficients for the linear profile case 
at 5 = 0. The values used for this normalization aic K„ = 44975 kN/m (256883 
Ib/in), C„ = 21.78 kN-s/m (124.4 Ib-s/in) and = 15821 kN/m (90364 Ib/in). 

For this study, the mid-point clearance of the quadratic profile is made 75% of 
(ci + Ce)/2, i.e., 0.75 times the mid-point clearance of a linear profile with «iTnila.r inlet 
and exit clearances. 

The plot for direct stiffness in Figure 5.8 shows the effect of a change in profile on 
the direct stififfiess. For the linear case, there is a complete loss of stifiness at around 
6 = 0.65. The stiffness for the quadratic profile is almost twice that of the 1inf»Rr 
profile. Also, it retains its stifiBaess over a much wider range than the linear profile. 
The difference in the other coefficients, shown in Figures 5.9—5.10, are relatively small. 
There is a drop of about 25% in leakage for the curved profile. 

Based on these results, the following conclusions may be drawn. 

1. Complete loss of sti&ess for linear profile at 6 = 0.65. 

2. The direct stiffiaess of curved profile is almost double that of linear profile. 

3. There is a 25% reduction in leakage for curved profile. 

This example illustrates the effect of seal profile on the dynamic coefficients. This 
example shows that and an arbitrary seal profile, other than a strai gh t or tapered. 
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could possibly be used 
characteristics. 


criterion in designing a seal for a set of optimum dyn ami r 


5.3 Case Study of a Distorted Seal of SSME-ATD-HPOTP 


Th. predicted cl«raoc. profUe of a. iotcratage of tha SSME-ATD-HPOTP t«r- 

bop^p i. .ho™ ia Figaae U. Tha diafoatad claarapca p^ffla is obtaiaad from a 

thaam«Ia.t.c &d.a alamant analysi. of tha turbopmnp. Tba claarm.ca. axa obtainad 

af a,mdi.tant axial plana, along tha langth of tha .aal with 6S daaranca. at aach 

Plana. Tha daaaanca. along tha cfrcntnfatanca am loaatad, tonghly, at agnal angnlat 
displacements. 

Tba ganaral procadnra amployad at NASA/MSFC with tha.a di.tortad profila. 
.. to compnta tha avataga inlat and axit alaaranca. and nmd than. a. inlat and axit 
anca. of a taparad Mai. In thi. .tndy, rotordynamic coeffidant. of tha di.tortad 
.«1 ata compatad with tho.a compntad naing avaxaga daatanca. at inlat and ont- 
l.t. Tha .aal gaomatry and oparating condition, at fell powar lavd (FPL) am givan 
m Appendix G. Tha daaranca function for thi. Mai i. approximatad by fitting tha 
craranca data with bi-cnbic .plinas. Thi. 2-D curve fitting enable, tha numerical 
cumputatton of daaranca. c(r,d), and gradient, g, fe at any given grid location 
According to tha manufactumr’. .pacification., tha .ida-load on thi. Mai act. 

.t a conatant angle of 290«. Tha dynamic coeffidant. for thi. vmiabla profila .aal are 
computed a. a function of .ida-load acting at thi. angle. Tha concept of axtamal load 
based dynamic coefficients is discussed in section 2.10. 

Kgura 6.12 .how. tha tha relation between Mai forc« m.d accantridty. At zero 
oad, dratortad profila .how. an accantridty. No load operation require, tha .aal to be 
chghtly offi-cantarad due to tha uneven di.tributmn of fluid praMura in tha di.tortad 



Normallzsd Crot* Coupled Stiffness, Kxy. Kyx 
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Figure 5.8 Normalized Direct Stiffness for Elliptical Seal (Curved) 



Figure 5.9 Normalized Cross Coupled Stiffness tor Elliptical Seal (Curved) 
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fo. the av„.,e and tha distorted prtdUe, is sH„„a in Fi^e 5.13. 

here .s smaU mcrease ip leaiage the distorted seal. 

Figure 5.14 shows the variation of direct stifeess K K 
_ , , , , stumcss, Kyy as 6. function of 

uxtumal load for both average and distorted proaes The cr™,, ,n 

sh™ • r- Proiues. ibe cross coupled stiffness 

Ml Figure 5.15 clearly shows the difference between the gg 

uua..is and astorted seal anapsis There ' “ 

at high loads. “ 

5.4 Directiona Dependence of Djmamic CoeBcients 

. in the case of a plain jonma bearing, the dpnandc coeffirtents of an 
aaa ar s are computed m a minimum am tadtness coordinate system (x' p-) 
.^hown a Figure 5.1f. The sea repres«ted in this hgure is of a rtrcuW Is 
2 -im to a plain jouma hearing. In this dgure, (x.p) represent, the globa 
coor^ate system hned to the stator and is the coordinate system normally used 
ro^rdynamic simulations. Therefore, the dynamic coeffidents irrespective of the 
coordmate system in w^ch they are computed shoad be transformed into ta. lined 
coordmate system befom they can be used in simuiaion stuaes. 

tridt • • A KMuputed at eccentridties greater than sero, th eccen- 

y IS vane dong the x'-anis of the minimum film thickness system and the 
compu e coeffiaent. am then transformed into the fixed coordinate system using a 
Wo^^on. For a sea with a circular cross section, dynamic coefiidents com- 
P-d a eccentridties greater than aero need to be transformed and those computed 
•t xero eccentridty need no transformation. 

In Figure 5.17, the angle of rotation between lb. ™ e. 

... .. “-umum film thickness system 

Mid the fixed coordinate system » A, mnsi • i ^ 

system is <f> and is also known as the eccentricity angle as 


Leakage. Q (kg/a) Eccentricity 


101 



Figure 5.12 Load vs. Eccentricity for Distorted Seal, Unit 3-01 



Figure 5.13 Leakage for Distorted Seal, Unit 3—01 





Crot* Coupled Stiffness, Kxy, Kyx (MN/m) 
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Figure 5.14 Direct Stiffness for Distorted Seal, Unit 3-01 



Figure 5.15 Cross Coupled Stiffness for Distorted Seal, Unit 3-01 
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Figure 5.16 Direct Damping for Distorted Seal, Unit 3-01 

it th, «c«.tricity vector with respect to the Sxed coordirutte systeot. The 

points O end C7 refer to the center of seal and rotor respectively and the line passing 
through these two points is the line of centers. 

The transformation of dynamic coefficients computed in a minimum film thick- 
ness system mto fixed coordinate system is explained below. 

Let (<?] be the transformation matrix between the minimum film thickness system 
specified by (x',p') coordinate system and and the global coordinate system repr^ 
seated by (i,p). The angle of rotation between these two coordinate systems is 
which m aUo the eccentricity angle as shown in Figure 5.18. Given a dynamic coef- 

fident in (x',p') system, it is required to Worm it into (z,y) coordinate system, 
which is the coordinate system nonnaUy used for dn.ul.n„.. 


0: center of eeal 
O’: center of rotor 


0—X’: line of centers 
9 '. *0000100117 angle 


Figure 5.17 Definition of Angle of Line 


of Centers and Eccentricity Angl( 



■jstem 


Figure 5.18 Rotor and Minimum F ilm 


Thickness Coordinate Systems 
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The transformation matrix [Q] and its inverse [Q] ^ is given below. 

cosS sin<f> 

[Q] = 

—sin<f> co8<f> 


(5.14) 


(5.15) 


cos<f> —sin<f> . 

[(}]-> = (5.15) 

sin<f> coB<j> 

m ■ 

Let the set of 12 dynamic coefficients at any given eccentricity e in the minim um 
filTTi thickness coordinate system (x\ y^) be specified by the stifihess matrix \K'\^ 
damping matrix [C"] and the inertia matrix [M'\. 


[r(e,0)] = 


[C"(e,0)] = 


[M'(e,0)l = 


K' 

XX 



KLy_ 

CL, 



CLy^ 

AC 



AC 


(5.16) 


(5.17) 


(5.18) 


Let the set of 12 dynamic coefficients at the same eccentricity e as above, but in 
the global coordinate system (x,y) be given by by the stififfiess matrix [iif], damping 
matrix [C] and the inertia matrix [Afj. 


[KM)] = 


K„ Ky 

— ky, Kyy 
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[C(e,«i)l = 



Cry 


Cyy 



-TUy, 

Myy 


( 5 . 20 ) 


( 5 . 21 ) 


[Af(e,^)l = 

Tie two Kts of dynamic coefficient, are related by the following tranrformation. 

me,4>)] = |(?l-Iiir'(e,0))[(J) (5.22) 

[C(e,^)l = |(?)-(C7'(e,0)I((?J (5.23) 

(A/(e,^)] = (<?)-‘(Jl/'(e,0))[(?] (5.24) 

For a «al with a circniar cross section, the coefficients compnted in the 

film thickness system will be the same irrespective of the eccentricity angle 4 ,. In 
other words, coefficients computed at two different orientations ^ and at * in the 
minimum fflm thickness system wffl he the same. Typically, the (x',p') is aUgned 

(*. p), l.e, ^ = 0, when these coefficients ate computed. The results given in seal 
literature are usually computed in this fashion. 

However, this procedure is no not valid with seals of non-circular cross sections, 
i.e., seals with circumferentially varying clearance functions. For these s«ds, the 
dynamic coefficients computed in the minimum film thickness vary with eccentricity 
angle m>d because of this the orientation of the (n',p') system with respect to (»,p) 
■s important. One way to handle this is to compute these coefficient directly in the 
fixed coordinate system for a given eccentricity angle, or compute the coefficients in 

the minimnm film thickness system and transform them using the transformations 
given in Eqs. (5.22-5.24), 



CroM Coupled Stlffneee, Kxy (MN/m) Direct 8tlffr<eee, Kxx (MN/m) 
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Figure 5.19 Variation of for Elliptical Seal 



Figure 5.20 Variation of for Elliptical Seal 
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Figure 5.21 Variation of C7«« for Elliptical Seal 

The following example illustrates the importance of orientation for arbitrary pro- 
file seals. Dynamic coefficients for a circular seal (ellipticity, ^ = 0) and an elliptical 
seal (elhpticity, 6 = 0.4), are computed using the two different approaches discussed 
earlier, for concentric position. Comparison plots for direct stiffiiess cross cou- 
pled stiffiiess and direct damping C„ computed using two approaches discussed 
above are shown in Figures 5.19-5.21- 

In these figures, Actual Cotfficitfits refer to coefficients computed directly in the 
fixed coordinate system (x^y) for a given eccentricity angle. The eccentricity an- 
gle is swept form 0 to 360® and these coefficients are computed at regular intervals. 
The Transformtd Coefficients are coefficients computed in a TniTiimtiTn film thick- 
ness system aligned with fixed coordinate system, and then transformed for a given 
eccentricity angle using the transformations in Eqs (5.22-5.24). 

For the 6^ = 0 case, i.e, for a circular seal, there is no difference between the 
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two approaches. The transformed coefficients coincide with actual coefficients at all 
orientations, as expected. The same procedure repeated for an elliptical seal (5 = 0.4), 
shows the differences between the two approaches. The transformed coefficients are 
entirely different from actual coefficients computed in the fixed coordinate system. 

Therefore, for seals ¥nth non*circular cross sections, either the coefficients are 
to be computed in the fixed coordinate system at a given eccentricity angle, or if a 
TTiiTiiTmim film thickness system is used, the actual eccentricity angle should be used 
to transform them into the fixed coordinate system. In other words, for these seals, 
dynamic coefficients should always be referred with respect to the orientation at which 
they are computed. 

This the primary reason for using a fixed coordinate system for the analysis of 
arbitrary profile seals, for example, the distorted seal case and the elliptical seal case. 
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CHAPTER VI 
RESULTS 

In this chapter, ciurent analysis is compared with both experimental and theo- 
retical results from literature. 

The following cases are studied in detail. 

1. Childs and Lmdsey (1993): This study presents experimental and theoretical 

results for high speed, short length, smooth, liquid annular seals with an axial 

taper. Experimental results are presented for both concentric and eccentric 

tests. Theoretical results for concentric case are based on Childs’ (1993) code 

MUDY, while similar results for eccentric tests are based on San Andres’ (1991) 

seal code HSEAL. The friction model is Moody’s and constant properties are 
assumed. 

2. Childs and Kim (1985): Theoretical and test results for a concentric seal based 
on Hirs’ friction model and constant properties. 

3. Sch«,^ and (1989): TWtical rasnlt. for a with a wavy (distorted) 

profo. in the axial direction. The «ction factor it ha«:d on Hire’ model and 
constant properties are used. 

4. Scharrer and Nelton (1990): Theoretical r«:nlt. for a parMy tapered annular 

>«d. Resnlts are for a concentric with Hirt’ friction model and conttant 
properties. 

5. Jenssen (1970): Experimental results (seal forces) for smooth annular seals as 
a function of eccentricity. 
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6. Kaaki and Kawakami (1984): Experimental results for long pump anTnilar seals. 
Theoretical predictions of Nelson and Nguyen are included. 

7. Falco et al. (1984): Experimental and theoretical results for plain «.TiTiiiljir seals. 
Theoretical predictions of Nelson and Nguyen are included. 

8. Allaire et al. (1976): Theoretical results based on short seal assumption and a 
Blassius type friction model. TheoreticiJ predictions of Nelson and Nguyen are 
included. 

9. San Andres et al. (1992): Theoreticsd resTilts for a cryogenic seal with Moody’s 
friction model, Isothermal flow with variable properties. 

10. San Andres et al. (1992): Theoretical results for a cryogenic seal with Moody’s 
friction model, Adiabatic flow with variable properties. 

6.1 Childs and Lindsey (1993) 

The results discussed in this section are based on the combined experimental and 
theoretical work of Childs and Lindsey (1993). A summary of this work is presented 
below. 

6.1.1 Work Summary 

This work presents theoretical and experimental results for water lubricated, short 
length, smooth, liquid a nn u l a r seals with an axial taper. Experiments are con- 
ducted with five different seal configurations at three different pressure differentials, 
1.38 Mpa, 2.41 Mpa, and 3.45 Mpa. The experiments are repeated at three speeds 
10200 rpm, 17400 rpm and 26400 rpm. 
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Table 6.1 Seal Geometry for Childs and Lindsey 


seal no. 

taper par. (q) 

c,(mm) 

Ce(mm) 

c.(mm) 

1, maximum divergent seal 

-0.29 

0.076 

0.137 

0.076 

2, slightly divergent seal 

-0.12 

0.076 

0.097 

0.076 

3, straight seal 

0.00 

0.076 

0.076 

0.076 

4, slightly convergent seal 

0.12 

0.097 

0.076 

0.076 

5, maximum convergent seal 

0.29 

0.137 

0.076 

0.076 


In this study, the experimental results are compared with the theoretical predic- 
tion of Childs’ (1993) computer code MUDY for concentric test results and with the 
predictions of San Andres’ (1991) seal code HSEAL for eccentric test results. San 
Andres analysis employs a hnite difference based solution scheme while Childs uses 
direct integration. The friction model used is Moody’s and constant fluid properties 
are assumed. 

The taper ratio of a seal is specified by the taper parameter q which is defined 


as, 



9 

and for, 


9 = 0 

straight seal geometry 

9 > 0 

convergent seal geometry 

9 < 0 

divergent seal geometry 


Cj-C^ 
Ci + Ct 


( 6 . 1 ) 


The five seal configurations used in the study are identified by their taper pa- 
rameter q as given in Table 6.1. 
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It is assumed in this study that the pre-swirl ratio for a short ATTnuUr seal is 
approximately equal to its whirl-frequency ratio. This ratio for each case is determined 
from the experimentally measured rotordynamic coefficients as, 


WFR = 

{Cgg + Cyy)u 


( 6 . 2 ) 


and are tabulated for all cases. This data is included in Appendix G. 

The seals are classified as smooth seals and the stator and rotor relative rough- 
ness is based on Moody friction model. The stator and rotor relative roughness, inlet 
loss coefficient, and exit pressure recovery coefficient are selected to match theoretical 
and experimental flow rates and the best set of values used in the theoretical predic- 


tions for all cases are given as, 
inlet loss coefficient, 0.1 

exit pressure recovery coefficient, 1.0 

stator rel. roughness (Moody, smooth) 0.001 

rotor rel. roughness (Moody, smooth) 0.001 


Experimental results are provided for concentric and eccentric seal tests. For 
the concentric position runs, all five seal configurations are used. For the eccentric 
position runs, results are provided only for the straight seal {q = 0) and the slightly 
convergent seal (q = 0.12). The experimental results include measured flow rates and 
rotordynamic coefficients. 

6.1.2 Comparative Study 

The experimental and theoretical data included in this study present an excellent 
opportumty to compare the results of the present work with the analyses of Childs’ 
and San Andres’ since one of the problems associated with studies involving combined 
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experimental/theoretical data is a lack of consistent input data for theoretical predic- 
tions and each researcher typicaUy chooses his or her own set of input data to match 
the experimental data. The objectives of this comparative study are two fold. One is 
to to compare the theoretical predictions of current analysis with experimental data 
and the other more important objective is to compare current analysis with Childs’ 
and San Andres’ analyses under si milar assumed input data. This comparative study 
IS significant in the sense that the results from the current analysis are being com- 
pared to the analyses of Childs and San Andres, who use solution procedures different 
from the current work, though all three analyses essentially use the same bulk flow 
governing equations, friction factor and boundary conditions. 

In the following comparisons, the theoretical predictions are repeated based on 
the current analysis using exactly the same input parameters, i.e., same pre-swirl, inlet 
loss coefficient, exit pressure recovery coefficient, stator and rotor relative roughness, 
density and viscosity as reported by Lindsey (1993). The predictions from current 
analysis are repeated for both nominal clearances and measured clearances. Measured 
clearances are clearances measured under running conditions and take into account 

(Lindsey, 1993) the rotor growth and change in nominal seal clearance during the 
operation. 

As will be evident from the comparative study to follow, the results based on 
nominal clearances (NCLR) are consistently closer to the experimental data than 
the results based on measured clearances (MCLR), suggesting possibly a need for 
refinement of the measurement system used for measuring these clearances. 

In the plots shown for this study, N refers to results based on nominal clear- 
ances, and M refers to results based on measured clearances. Nominal clearances 
based results are available only for current analysis. For comparisons between various 
analyses, i.e, current analysis, Childs and San Andres, the results based on measured 
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clearances are used. 

In the following study, results from current analysis are provided for all cases 
given below. 

1. 5 seal configurations, taper par. q: -0.29, -0.12, 0.00, 0.12, 0.29 

2. 3 pressure differentials, Ap: 1.38 Mpa, 2.41 Mpa, 3.45 Mpa 

3. 3 rotor speeds, N: 10200 rpm, 17400 rpm, 24600 rpm 

4. Concentric and eccentric tests. 

6.1.3 Leakage, Concentric Tests 

The leakage for concentric seal operation as a function of taper parameter q is 
given m Figure 6.1. As mentioned earlier, Lindsey uses Childs’ code MUDY based 
on Moody’s friction model, for theoretical predictions of concentric seal tests. These 
predictions from MUDY along with the predictions from the current analysis for both 
NCLR and MCLR are shown in this figure. As may be noted from this plot, there is 
a considerable difference between predictions based on MCLR and the experimental 
data. However, the leakage based on NCLR show very good agreement particularly 
for the convergent seal geometry where there is almost exact correlation with the 
experimental data for all speeds and pressure differentials. The muTiTmiTn deviation 
IS about 12% and it occurs for the maximum divergent case at 3.45 pressure differential 
and 24600 rpm case. This result is directly opposite to the results based on MCLR 
which are closest to experimental data for the maximum divergent case. In almost 
all cases, current analysis with MCLR predicts leakage which is slightly (about 10- 
15/Childs’ predictions. One of the possible reasons for current analysis based on 
NCLR being much closer to the convergent seals’ results compared to the divergent 



iMkao* (Nlcrt/mln) 


6 



taper paraiMlar (q) 


Figure 6.1 Leakage, Childs and Lindsey, Concentric Tests 
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seals may be due to flow separation that is likely to occur in divergent seal geometries 

(Scharrer and Nelson, 1990) and bulk flow model used is not equipped to deal with 
that type of flow. 


6.1.4 Dynamic Coefficients, Concentric Tests 


The rotor dynamic coefficients for the concentric tests are shown in Figure 6.2. 
for 10200 rpm case and in Figure 6.3 for the 17400 rpm case. The correlation between 
theoretical and experimental data for direct stiffiiess is, at best, average. However, 
the present analysis results based on NCLR show a much better correlation with 
experimental data than the MCLR based analysis, particularly for higher pressure 
differentials. Both current analysis (MCLR) and Childs agree weD with the experi- 
mental data for the highly convergent case. Also, the theoretical predictions generally 
follow the trend of the experimental results, i.e., increase in sti&ess with taper param- 
eter q. The maximum deviation for direct strffiiess occurs for the TnA-riTnnm divergent 
case and a possible reason is flow separation as mentioned earlier. 

Both MUDY and current analysis predict simUar damping and cross coupled 
stiflbess. The added mass is severely under-predicted by both analyses, and Lindsey 
(1993) points out that this big difference may be due to unaccoimted fluid inertia 
effects in the housing and piping system. In spite of the above reason, theoretical 
predictions typically under-predict mass coefficients. 

Similar trends are noted for the dyneonic coefficients of the 17400 rpm ex- 
cept for the damping where the difference between theoretical and experimental data 


mcrease. 



Direct Damping. Cxx (kN>t/m) OIrgct Stlffnata. Kxx (MM/m) 
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Figure 6.2 Dynamic Coefficients, Childs and Lindsey, Concentric (10200 rpm) 






OIrtcl Dampino. Cia (kN-t/m) OIract StlfVnaM. Kn (MH/m) 
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Figure 6.3 Dynamic Coefficients, Childs and Lindsey, Concentric (17400 rpm) 
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6.1.5 Leakage, Straight Seal, Eccentric Tests 

For eccentric tests, experimental data is available only for two seal configurations; 
a straight seal (9 = 0) and a slightly convergent seal (9 = 0.12). The code used for 
theoretical predictions for eccentric seal analysis is HSEAL, developed by San Andres 
(1991) and is based on a finite difference formulation. The plots shown in Figure 6.4 
correspond to the leakage of a straight seal operated at an eccentric position with 
eccentricity ratios varying from 0 to 0.5 and for three speeds. The prediction of the 
leakage rate by current analysis is typically 10-15% better than HSEAL predictions. 

However, the maximum difference between test data and current analysis based on 
NCLR is only about 10%. 

The large deviations predicted by HSEAL are inexplicable and 

6.1.6 Dynamic Coefficients, Straight Seal, Eccentric Tests 

The plots in Figure 6.5 refer to the dynamic coefficients for the straight seal op- 
erated at various eccentricities. Predictions by both the current analysis and HSEAL 
are similar and follow the trends of the experimental data. Good comparison for 
direct damping and cross coupled stiffiiess for both analyses. 

Similar trends are seen for the 17400 rpm case. 

6.1.7 Leakage, Slightly Convergent Seal, Eccentric 

The results in this section correspond to the leakage of the shghtly convergent seal 
{q = 0.12). As noted earlier, there is excellent correlation between current analysis 
(NCLR) and measured flow rates for the convergent geometry seals (Figure 6.1). 
This very good correUtion is repeated for the case of sKghtly convergent seal shown 
in Figure 6.7. Best comparison of flow rates for all seal configurations tested. These 
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Figure 6.4 Leakage, Childs and Lindsey, Straight Seal, Eccentric Tests 






Crott Couptod SmrnMt, luy (MNAn) Dlr«« SUffTMti, Kn (iM/m) 
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Figure 6.5 Dynamic Coefficients, Cliilds and Lindsey, Straight Seal, Eccentric (10200 
rpm) 
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Figure 6.6 Dynamic Coefficients, Childs and Lindsey, Straight Seal, Eccentric (17400 
rpm) 
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results based on nominal dearances point to the fact that possibly the measurement 
system employed to measure dearances during operation need to be refined. 

There is considerable difference of more than 20% in leakage predictions between 

current analysis (MCLR) and San Andres’ results, particularly for 17400 rpm and 
24600 rpm cases. 

6 . 1.8 Dynamic Coeffidents, Convergent Seal, Eccentric Tests 

The plots m Figure 6.8 and Figure 6.9 correspond to the dynamic coeffidents for the 
slightly convergent seal configuration at 10200 rpm and 17400 rpm respectively. Ex- 
cellent companson of flow rates (NCL) translate into better correlation with dynamic 
coeffiaents, particularly for direct stiffness. For current analysis based on MCLR 
and San Andres report similar dynamic coeffidents, with a slight variation in direct 
sti&ess for the higher pressure differential case. 

6.1.9 Condusions 

Lindsey makes the following condusions based on the above study. 

1. In general, results are consistent with theoretical predictions, except for flow 
rates. 

2 . Theory largdy under predicts flow rates. 

3. Flow rates comparison best for maximum divergent case 

4. Direct Stiffiiess increases with taper parameter q. 


5. Cross coupled stiffness predicted well by theory. 

6 . Damping decreases for 9 < 0 and 9 > 0 . 



LMk.O« (lll.fi/mln) (IH^Wmln) 
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Figure 6.7 Leakage, Childs and Lindsey, Convergent Seal, Eccentric Tests 
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Figure 6.8 Dynamic Coefficrents, Childs and Lindsey, Convergent Seal, Eccentric 
(10200 rpm) 
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Figure 6.9 Dynamic Coefficients, Childs and Lindsey. Convergent Seal, Eccentric 
(17400 rpm) 
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7. Damping increases with eccentricity ratio. 

8. For eccentric seals, the dynamic coefficients remain relatively constant upto an 
eccentricity ratio of 0.5. 

Current analysis predictions based on nominal clearances for flow rates are much 
closer to measured flow rates than measured clearances based analyses. For conver- 
gent seal geometry, the flow rates almost match exactly. When measure clearances 
are used, current analysis shows a slightly better agreement (10-15%) than either 
MUDY or HSEAL. Comparison of all analyses in the case of direct stifhiess is, at 
best, average. Current analysis based on NCLR gives a better correlation than the 
other analyses based on MCLR. 

6.2 Childs and Kim (1985), Hirs Model 

Childs and Kim [1985] presented an analytical and experimental study for rotordy- 
namic coefficients of turbulent annular seals with different directionally-homogeneous 
surface roughness treatments for rotor and stator surfaces. The friction model is 
based on Hirs’ model and constant properties are assumed and the analysis is for a 
concentric seal. The seal code based on the the analytical part of this study had been 
the mainstay of seal analysis work at NASA/MSFC for many years. MSFC provided 
this author with data for some test cases which were then compared with the results 
from the present analysis for Hirs’ friction model. Results from one of the test cases 
is given in Table 6.2. The seal data for this example is given in Appendix G. 

The results from current analysis for Hirs’ friction model match well with the 
results of Childs and Kim (1985). 
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Table 6.2 Cbilds and Kim Check Case, Hirs’ Model 


Data 

Ckilds/Kim 

Current Analysis 


88.14 MN/m 

88.14 MN/m 

kmy 

11.11 MN/m 

11.11 MN/m 

c„ 

16.15 kN-s/m 

16.14 kN-s/m 

Cry 

0.419 kN-s/m 

0.419 kN-s/m 


0.3215 kg 

0.3214 kg 

^ry 

-0.0024 kg 

0.0005 kg 

Q 

7.959 kg/s 

7.959 kg/s 


6.3 Scharrer and Nunez (1989) 

The effect of seal distortions on rotordynamic coefficients was first considered by 
Sharrer and Nunez (1989). They reported that a 2-D, axisymmetric, finite element 
analysis which considered the intemsil pressure distribution, and the boundary condi- 
tions due to assembly and operating interferences produced a clearance profile which 
was wavy and different from the nominal design tapered profile. 

This distorted seal profile in the axial direction was fitted with a clearance func- 
tion in the form of a polynomial as, 

h{z) = oi -f Q 2 Z -f osz’ -I- 04 Z* + aiz* (6.3) 

where the coefficients Oi, 02 , • • • etc., are coefficients chosen to fit the distorted axial 
profile. 

They adapted the analysis of a plain seal to the case of a wavy profile seal. They 
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Table 6.3 Scharrer and Nunez, Rx>ugh Wavy Seal Case 


Data 

Scharrer /Nunez 

Current Analysis 


46.35 MN/m 

57.30 MN/m 


52.02 MN/m 

51.57 MN/m 

c„ 

33.42 kN-s/m 

33.65 kN-s/m 


1.49 kN-s/m 

1.43 kN-s/m 


0.753 kg 

0.772 kg 


0.026 kg 

0.015 kg 

Q 

0.471 kg/s 

0.462 kg/s 


reported a marked change in the computed rotord 3 rnamic coefficients due to a change 
in the seal profile. These changes include, a) loss in direct stiffiiess b) increase in 
cross coupled stiffiiess c) increase in damping. The results for the case of rough wavy 

seal is given in Table 6.3. The direct stiffiiess between two analyses differ by about 

20 %. 


6.4 Scharrer and Nelson (1990) 

Scharrer and Nelson (1990) conducted a theoretical study of an MtinlAr seal with a 
partially tapered clearance. In this study, they investigated the axial distortion prob- 
lem. They tried to correct the predicted distortions by machining out the undesirable 
distortions at the design stage itself. The model they used to accomplish this is a 
seal with a taper on part of length of the seal. Using this model, they conducted a 
parametric study of various performance characteristics as a function of taper length 
to total length ratio (T/L). Based on this study they recommended optimum ratio 
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Figure 6.10 Scharrer and Nelson, A Partially Tapered Annular Seal 

of T/L for best performance of these partially tapered seals from a rotordynamic 
analysis point of view. 

They developed the analysis based on Hirs’ turbulent lubrication equations (Hirs’ 
friction model). In this analysis, the seal is assumed to have a taper only over a portion 
of the length of the seal and the rest of the seal is treated as a straight seal. 

Figure 6.7 shows the details of a partially tapered seal. L is the total length 
and T is the taper length and and c, are inlet and exit clearances. The taper 
length/total length ratio is varied by varying the parameter, q given as, 

9 = TjL (6.4) 

The seal is a completely straight seal for g = 0 i.e., c, = c, and a fully tapered seal for 
9 ~ ^*6* This ratio is varied from 0 to 1 smd its effect on leakage and rotord}maniic 
coefficients is studied. 

They analyzed two different seal configurations based on stator and rotor relative 
roughness. 

1. Taper Smooth (varying g) 

2. Taper Rough (varying g) 
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Table 6.4 Scbarrer and Nelson, Partially Tapered Seals (Smooth), Taper=1.0 


Data 

Scharrer /Nelson 

Current Analysis 

K„ 

242.0 MN/m 

215.0 MN/m 

Ky 

65.5 MN/m 

62.6 MN/m 

c„ 

26.0 kN-s/m 

24.8 kN-s/m 


0.66 kg 

0.66 kg 


For both the above cases, the inlet clearance to cadt clearance ratio ^ is main- 
tained at 3. A portion of results from that study along with current analysis results 
are presented in Tables 6.4-6.Q. 

6.4.1 Smooth Seals 

For the case of smooth seals, results are compared for three taper ratios; 

1. taper ratio, g = 1, a fully tapered seal. 

2. taper ratio, ^ = 0.4, a partially tapered seal. 

3. taper ratio, q = 0.0, a fuUy straight seal. 

Results given in Tables 6.4— 6.6 show about 20% difference in and slightly 
smaller deviations for 1^, for all cases. 

6.4.2 Rough Seals 


Similar results are reproduced for the case of rough seal for three different taper 
ratios. The results are shown in Tables 6.7— 6.9. Again, as in the case of smooth seals. 
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Table 6.5 Scharrer and Nelson, Partially Tapered Seals (Smooth), Taper=0.4 


Data 

Schamr/N elson 

Current Analysis 

K„ 

220.0 MN/m 

179.7 MN/m 

Ky 

82.0 MN/m 

75.26 MN/m 

c„ 

32.0 kN-s/m 

29.45 kN-s/m 

M„ 

0.55 kg 

0.55 kg 


Table 6.6 Scharrer and Nelson, Partially Tapered Seals (Smooth), Taper=0.0 


Data 

Scharrer/Ndson 

Current Analysis 

K„ 

152.0 MN/m 

159.0 MN/m 

/Sey 

102.5 MN/m 

105.2 MN/m 


39.2 kN-s/m 

40.1 kN-s/m 


1.02 kg 

1.02 kg 
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Table 6.7 Schatxer and Nelson, Partially Tapered Seals (Rough), Taper=1.0 


Data 

Scharrer /N elson 

Current Analysis 


222.0 MN/m 

194.0 MN/m 

kfy 

61.0 MN/m 

57.9 MN/m 


25.8 kN-s/m 

24.3 kN-s/m 

Af*x 

0.66 kg 

0.66 kg 


difference of about 20% in direct stiffness. Both analyses use same governing 
equations with Hirs’ friction model. 

To explain this discrepancy, the results from this study are compared to Childs 
and Kim (1985) as all three analyses use the same governing equations based on 
Hirs friction model and differing only in the solution procedure adopted. The results 
of this comparison are shown in Tables 6.10-6.11. These results show Scharrer and 
Nelson’s analysis differs from Childs’ and current analysis consistently. It is likely 
that Scharrer was using the same analysis that he and Nunez (1989) used for the 
wavy profile seal analysis where a similar discrepancy was also noted. The seal data 
for this study is included in Appendix G. 

6.5 Jenssen (1970) 

Jenssen (1970) investigated experimentally the load bearing capacity of smooth liquid 
annular seals at various eccentricities. The test data was collected for three pressure 
differentials, 0.344 Mpa, 1.034 Mpa, 1.724 Mpa and at three different speeds 3000 rpm, 
5000 rpm and 7000 rpm. The seals' used in the experiment are long seals (L/D=1.025) 
with water as the working fluid. He also presents theoretical predictions based on 
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Table 6.8 Scbarrer and Nelson, Partially Tapered Seals (Rough), Taper=0.4 


Data 

Schamr/N tlson 

Current Analysis 


190.0 MN/m 

170.2 MN/m 


68.0 MN/m 

70.0 MN/m 


29.0 kN-s/m 

29.0 kN-s/m 


0.67 kg 

0.59 kg 


Table 6.9 Scbarrer and Nelson, Partially Tapered Seals (Rough), Taper=0.0 


Data 

Scharrer/Nelson 

Current Analysis 


118.0 MN/m 

121.0 MN/m 

Kv 

87.5 MN/m 

89.3 MN/m 

c„ 

38.2 kN-s/m 

38.6 kN-s/m 


1.08 kg 

1.09 kg 
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Table 6.10 Scharrer and Nelson, Comparison with Childs and Kim, Tapcr=1.0, Rough 


Data 

Scharrer/Nelson 

Childs/Kim 

Current Analysis 


242.0 MN/m 

215.7 MN/m 

214.8 MN/m 


65.52 MN/m 

62.74 MN/m 

62.69 MN/m 

c„ 

26.0 kN-s/m 

24.8 kN-s/m 

24.8 kN-s/m 


- 

3.16 kN-s/m 

3.16 kN-s/m 


0.66 kg 

0.66 kg 

0.66 kg 

TTlgy 

- 

0.0027 kg 

-0.0041 kg 

Q 

- 

9.568 kg/s 

9.551 kg/s 


Table 6.11 Scharrer and Nelson, Comparison with Childs and Kim, Taper=1.0, Rough 


Data 

Scharrer /Nelson 

Childs /Kim 

Current Analysis 

K„ 

222.0 MN/m 

194.6 MN/m 

193.8 MN/m 


61.0 MN/m 

58.1 MN/m 

57.9 MN/m 


25.80 kN-s/m 

24.36 kN-s/m 

24.30 kN-s/m 


- 

2.72 kN-s/m 

2.72 kN-s/m 


0.66 kg 

0.66 kg 

0.66 kg 


- 

-0.0042 kg 

-0.0074 kg 

Q 

- 

8.356 kg/s 

8.346 kg/s 
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Figure 6.11 Seal Force at 3000 rpm, Jenssen 

shoH stal assumption which are not shown here. The seal input data for this case is 

taken from Nguyen (1988) and is given in the Appendix G. Results from the current 

« 

analysis along with the experimental data are given Figures 6.11-6.13, for the three 
rotor speeds. 

The steady state seal forces are plotted as a function of eccentricity ratio and 
the predictions from the current analysis agree well with the experimental data for 
all pressure differentials and at all speeds. 

6.6 Kanki and Kawakami (1984) 

Kanki and Kawakami (1984) investigated the dynamic bearing effects of long pump 
annular seals as a function of eccentricity. Nguyen (1988), reported convergence 
problems at eccentricity ratios of above 0.4 with the original method as shown in the 




S«al Force (N) S„, 
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Eccentricity Ratio 


Figure 6.12 Seal Force at 5000 rpm, Jenssen 



Eccentricity Ratio 


Figure 6.13 Seal Force at 7000 rpm, Jenssen 
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plots Figures 6.14-6.18. Beyond 0.4 eccentricity ratio, Nguyen’s method has problems 
converging and he attributes this convergence problem to to the onset of laminar flow 
and inability of the analysis to deal with negative direct sti&iess. However, with*the 
current analysis such a problem is not encountered. The test seal is a long smooth 
seal (L/D=1.0) and the working fluid is water. 

The results for this case along with Nguyen’s results are shown in Figures 6.14- 
6.18. It is true that the solution procedure has problems dealing with negative direct 
stifeess. However, this condition is rarely encountered in practical seal design and 
thus need not of major concern as far as the efficacy of this analysis is concerned. 

Figure 6.14 shows the steady state seal force as a function of eccentricity. Nguyen 
(1988) results are only upto 0.4 due to numerical problems while the current analysis 
gives reasonably good results upto 0.65. 

Nguyen s analysis predicts direct stiffiiess to decrease with eccentricity as 
shown in Figure 6.15. Current analysis predicts this direct stifihess in line with the 
test data, i.e., mcrease with eccentricity for Kyy and decrease with eccentricity for AT„. 
Also, Nguyen’s analysis runs into convergence problems around 0.4 eccentricity ratio. 
Current analysis encounters no such problems. There is reasonably good comparison 
for cross coupled stifeess, kyg, between current analysis and test data. 

The main feature of the results presented in this case is the better convergence 
properties of the current analysis based on cubic splines compared to the original 
approach of Nelson and Nguyen based on Fast Fourier Transforms (FFT). Also, the 
results from current analysis agree with experimental data better than Nelson and 
Nguyen’s analysis. 



Direct Stiffness, Kxx, Kyy (MN/m) 
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Figure 6.14 Seal Force, F, for Kanki and Kawakami 



Figure 6.15 Direct Stiffness, ii:„, Kyy for Kanki and Kawakami 








Direct Damping, Cxx. Cyy (kN-e/m) Crosa Coupled Stlffneta, Kxy, Kyx (MN/m) 
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Figure 6.16 Cross Coupled Stiffness, k^y, for Kanki and Kawakami 



Figure 6.17 Direct Damping, Cxi, for Kanki and Kawakami 
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Eeoantrlolty Ratio 

Figure 6.18 Cross Coupled Damping, c^y, Cy* for Kanki and Kawakami 
6.7 Falco ti al. (1984) 

In this combined experimental tind theoretical work, Falco et al. investigated the 
effect of eccentricity on the dynamic coefficients of an annular seal. Their analytical 
work was based on a iimte element model and they compared the analytical results 
with experimental test data. They compared their theoretical predictions with the 
various methods in use at that time and concluded that their finite element based 
analysis provided the best comparison with experimental data. Subsequently, Nguyen 
(1988) showed that the predictions &om his analysis were in better agreement than 
Falco’s theoretical results. In plots shown in Figures 6.22-6.26, Falco’s experimental 
results along with Nguyen’s predictions and current analysis results are given. The 
input data used is from Nguyen (1988) and is given in Appendix G. The results show 
good comparison with the original Nguyen’s approach. 



Cro«8 Coupled Stiffness, K*y. Kyx (MN/m) 
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Figure 6.19 Direct Stiffness K„, Kyy, for Falco et al. 



Figure 6.20 Cross Coupled Stiffness k^y, ky^ for Falco tt al. 







Direct Damping, Cxx. Cyy (kN-a/m) 
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Figure 6.21 Direct Damping, Cx*, for Falco et al. 



Figure 6.22 Direct Inertia, Mxx, for Falco et al 
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Figure 6.23 Leakage, Q, for Falco ti ol. 

6.8 Allaire et al. (1976) 

In this work, Allaire et al. investigated the effects of large eccentricity on the dy- 
namic coefficients of interstage seal of Space Shuttle Main Engine ffigh Pressure Fuel 
Turbopump (SSME-HPFTP). They used a solution approach using the short seal 
assumption (Couette flow) and a Blassius-type turbulent friction factor model. The 
results from current analysis along with Nelson and Nguyen predictions are given in 

Figures (6.19-6.21). The seal input data is taken from Nguyen (1988) and is included 
in the Appendix G. 


6.9 Comparison of Variable Properties Model with Constant Properties Model 

Results from a distorted seal analysis have been discussed in Chapter V. This exercise 
is repeated for the variable properties model developed in Chapter III. Comparisons 




Cross Coupled Stiffness, Kxy, Kyx (MN/m) Direct Stiffness, Kxx. Kyy (MN/m) 
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Figure 6.24 Direct Stiffness, iif„, for Allaire et al. 



Figure 6.25 Cross Coupled Stiffness, for Allaire ef al 
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Figure 6.26 Direct Damping, Cyy for Allaire ti al. 

made between the two modeb are shown in Figures 6.22-6.27. 

The main difference between these two models is a reduction in direct stiShess for 
the variable properties case and it difference with eccentricity as shown in Figure 6.22. 
This may be explained by the inclusion of compressibility into the analysis, in effect 
making the spring softer. There is negligible difference in other coefficients. 

Leakage, shown in Figure 6.27, is expectedly smaller for the variable properties 
model due to a decrease in density. 


6.10 San Andres et al. (1992), Isothermal Case 


San Andres ei al. (1992) presented theoretical results for a straight seal with thermal 
effects and variable fluid properties. Two cases are considered. 


1. Isothermal flow 



Cross Coupled Stiffness, Kxy, Kyx (MN/m) Direct Stiffness. Kxx, Kyy (MN/m) 
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Figure 6.27 Direct Stiffness K^, for Seal Unit 3-01 
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Figure 6.28 Cross Coupled Stiffness ik^,, for Seal Unit 3-01 
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Figure 6.29 Direct Damping, C^x, C^y for Seal Unit 3-01 



Figure 6.30 Direct Inertia, Afxx, Myy for Seal Unit 3-01 



Leakage (Kg/s) Seal Forca (N) 
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Figure 6.31 Seal Force, F, for Seal Unit 3-01 
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Figure 6.32 Leakage, Q, for Seal Unit 3^01 
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2. Adiabatic flow 

Tbe case of isotbennal flow (constant temperature) is the same as the variable 
properties model discussed in Chapter HI, i.e., fluid properties are assumed to be a 
function of local pressure and a mean temperature. The results for the isothermal 
case from current work are compared with San Andres’ in Figures 6.28-6.33. Direct 
stifihess, shown in Figure 6.28, difiers by about 20% and this difference is maintained 
at all eccentricities. Agreement for cross coupled stiflness, in Figure 6.29, is slightly 
better. 

6.11 San Andres et ai. (1992), Adiabatic Case 

The same case considered in the previous example is repeated with the thermal effects 
model of Chapter IV. The results are shown in Figures (6.34-6.38). are for the 
adiabatic case (Q, = 0, no heat transfer). The results are for concentric case and do 
not include perturbations in temperature as expired in Chapter IV. 

The temperature rise across the seal is shown in Figure 6.34. San Andres points 
out that if this rise in temperature is big enough, the fluid may enter a two-phase 
region seriously affecting the performance of the turbomachine. 

The frictional torque and flow rates shown in Figures 6.35,6.36 roughly match. 
However, cross coupled stiffness, shown in Figure 6.37 is off by about 25%. Damping, 
both C„ and agree well. 

6.12 Comparison of Current Analysis with Other Methods 

The following general conclusions may be drawn based on the check cases discussed. 
Comparison with Nelson and Nguyen: 

Expectedly, all check cases with current analysis compare well with their theoretical 



Croti Coupled StIHnost, Kxy. Kyx (MN/m) 
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Figure 6.33 Direct Stiffness Kxxi Kyy, for San Andres, Isothermal 



Figure 6.34 Cross Coupled Stiffness k^y, ky^, for San Andres, Isothermal 




Direct Damping, Cxx. Cyy (N-e/m) 
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Figure 6.35 Direct Damping, Cx*, Cyy for San Andres, IsothermaJ 



Eccentricity Ratio 


Figure 6.36 Direct Inertia, Mxx^ San Andres, Isothermal 




Ltakage (Kg/s) Seal Forca (N) 



Figure 6.37 Seal Force, F, for San Andres, Isothermal 



Eccantriolty Ratio 


Figure 6.38 Leakage, Q, for San Andres, Isothermal 
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Figure 6.39 Temperature Rise, San Andres, Adiabatic 



Figure 6.40 Frictional Torque, San Andres, Adiabatic 




StIMnett, Kxx. Kxy (MN/m) Leakage, (kg/e) 


156 


a.o - 
2.8 - 




■ 






2.6 - 


■ 

B 



1 





2.4 - 


MB 







2.2 - 










2.0 - 



4 






1.6 - 





'"i 




1.6 - 









1.4 - 









1.2 - 









1.0 - 










' I I I r 

10000 16000 20000 28000 30000 36000 40000 

Rotor Spaed (rpm) 


Figure 6.41 Leakage, San Andres, Adiabatic 



Figure 6.42 Stiffness, San Andres, Adiabatic 
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Figure 6.43 Damping, San Andres, Adiabatic 

predictions. In some comparisons with the experimental data, the comparison with 
current analysis is better than their results. For check cases where their method had 
convergence problems and failed, current analysis gives good results. 

Comparison with Childs: 

The current analysis compares very well with Childs’ Hirs’ friction model based anal- 
ysis (1985). For this friction model, current analysis matches Childs’ analysis almost 
exactly. Similar comparison exists between curren analysis and his more recent work 
based on Moody’s model. 

Comparison with Scharrer and Nelson: 

There is some deviation between their analysis based on Hirs’ friction model and the 
current analysis. A three-way comparative study between Childs’ analysis, Schar- 
rer and Nelson’ work and the current analysis shows discrepancy in their results. 
While Childs and current analysis agree well consistently, their results for stiffness 
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coefficients are off by about 20%. 

Comparison with San Andres: 

There is some deviation between the current analysis and San Andres’ Moody friction 
factor based eccentric analysis, both for constant properties and variable properties 
cases, particularly for stiffness coefficients. It appears as though these deviations vary 
from case to case. For example, for Childs and Lindsey (1993) experimental results, 
the difference in flow rates for convergent seals is considerably large. Same for direct 
sti&ess at high pressure differentials. However, for other cases the differences are not 
of that order. It is difficult to speculate on the reasons for these deviations as both 
analyses use entirely different solution procedures. 

However, for the current analysis an equivalence will be established between the 
d}mamic coefficients based transient motion and the same motion based on original 
governing equations. For the second approach there is no first order solution involved. 
If these two approaches match consistently, it establishes, in the minumum, that the 
is no error in the linearized coefficients obt^ed from the dynamic analysis. 



159 


CHAPTER VII 
TRANSIENT ANALYSIS 

In simnlating the dynamics of a rotor system with fluid bearings such as journal, 
tilt-pad bearings etc., or an a nnular seal in the present study, the dynamic efiects 
of seals for a small motion of the rotor about an equilibrium position are usually 
modeled using a lintariztd forct-moiion model Mmilar to the one shown in Eq. (7.1). 

In this equation, (^i,5y) are the displacements, (^x, jy) are the velocities and 
{Sx, Sj}) are the accelerations in the JC and Y directions respectively, relative to a 
static operating point (x,y). The fluid force terms AF, and AF^ are the incremental 
or perturbed fluid forces for a smaU motion of the rotor shaft about (x,y). These 
force components, in general, vary as a function of rotor displacement, translational 
velocity and acceleration and are linear only for small orbital motion. 

In this model, ky^ are the linearized sti&ess coefficients, 

c^, Cyt are the linearized damping coefficients and Af„, TTiy, are the lin- 

earized added mass or inertia coefficients at the static operating point or eccentricity 
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(7.1) 


In the linearized model, the terms [K„8z] and [KyySy] account for the incremen- 
tal fluid reaction forces of the s^ due to a small displacement of the rotor (5x,^y). 
The term is the cross coupled force in the X direction due to a displacement 6y 
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in the V direction. Similarly, [ky^Sx] is the cross coupled force in the Y direction due 
to a displacement Sx in the X direction. The terms \C^8x] and \Cyy8y] represent 
the incremental damping forces due to a small velocity change (tfz,5y). Similarly, 
\M„8x] and [Myy8y] are the incremental fluid inertia forces due to a small change in 
acceleration (8x,8y). For a concentric seal, K„ = - ky^ etc., reducing the 

number of coefficients from twelve to six. Typically, for an annular seal, the important 
coefficients are direct stiffness, cross coupled stiffness, direct damping and direct or 

added mass. The contributions of other terms are negligible in most cases compared 
to these terms. 

These twelve linearized coefficients are, in general, nonlinear functions of the 
static operating point (z,y). The variations of direct stiffiiess K„, direct damping 
C„ and cross coupled stiifriess k^ for various rotor operating positions for seal unit 3- 
02, an experimental seal under design at NASA/MSFC, are shown in Figures 7.2-7 A. 
These curves are obtained by the dividing the circumference of the seal into a number 
of segments as shown in Figure 7.1 and computing the coefficients as a function of 
eccentricity along each of the radii. For this seal, the coefficients remain constant 
upto an eccentricity ratio of 0.4. Beyond this limit, the coefficients start varying and 
this variation becomes much more pronounced as the eccentricity ratio exceeds 0.6. 
These curves are typical of a tapered seal and similar curves can be obtained for a 
straight seal. 

In practice, usually a single set of dynamic coefficients computed at centered 
position is used to model the dynamic behavior of the seal i.e., to compute seal 
farces in rotordynamic simulations. This is based on the experimental and theoretical 
observations (Childs 1993) that generally there is little change in dynamic coefficients 
upto an eccentricity ratio of 0.4-0.5. In other words, for simulations involving motion 
with in this range the dynamic coefficients computed at zero eccentricity should be 
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Figure 7.1 Circumferential Grid for Seal Coefficient Mapping 

adequate to model the seal behavior. For example, at NASA/MSFC for the SSME 
turbopump simulations, dynamic coefficients used to model the interstage seal are 
computed at zero eccentricity. 

It is assumed that this set of coefficients computed for a concentric seal would 
reliably predict the dynamic behavior of the seal over its entire range of operation, 
which may include motion with large eccentricities. This fact of large eccentric motion 
has been confirmed by the presence of destructive rubs in the SSME turbopump 
interstage seals. 

This method of modeling a seal using a single set of coefficients is valid only if 
the dynamic coefficients remain invariant in the clearance space. For example, for 
the seal unit 3-02 (Figs. 7.2-7.4) the coefficients remain relatively constant as long as 
the operating point falls within a circle of radius of about 0.4 eccentriaty ratio. As 
this limit is exceeded, the coefficients start varying and the variations are more rapid 
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at higher ecccntridties. 

In general, the above seal model (Eq. 7.1) consisting of 12 coefficients accurately 
approximates the dynamic behavior of the seal subject to a few limitations given 
below. 

1. The model is vahd only for a small motion in the immediate neighborhood of 
the static operating point at which the set of dynamic coefficients are computed, 
typically upto 0.4 eccentricity ratio. This is the basic assumption on which the 
linearized coefficients of the model are derived. 

2. The dynamic coefficients derived at a given static operating position may not 
be accurate when used at a different operating point. 

3. Even though these dynamic coefficients can be computed at various eccentrici- 
ties, in general it is not possible to decide which set of coefficients to use when 
the rotor is moving around in the clearance space such as in a transient motion. 


7.1 Objectives 

The main objective of this work is to study the effect of large rotor displacements 
of SSME-ATD-HPOTP turbopump on the d}mamic6 of the annular seal and the 
resulting transient motion. 

For the purpose of this study, large eccentric motion is classified into two types 
as illustrated in Figures 7.5-7.6. Figure 7.5 shows the time-displacement curve for 
the center of a rotor executing a steady state motion with a large amplitude. The 
amplitude is of the order of radial clearance (in this case about 0.017 mm) and hence 
naay be considered as a motion with large displacement. 

The motion represented in Figure 7.6 is of the second type where the rotor is 
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Figure 7.2 Vanation of for unit 3—02 as a function of rotor position 
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displaced to a large static eccentric position (c > 0.4) while executing motion gimilAr 
to type 1, The following study is not limited to these two types of motion and these 
arc used only for the purpose of illustration. The analysis to be developed is valid for 
any type of general motion. 

Results from the investigation of the first case should help in establishing the 
limits of accuracy of the linear force-motion model at a given static operating point. 
The basic xmderlying assumption of this model is that it is valid only for a small 
motion, typically for e > 0.4, about the operating point. The exact limits of this 
small motion are undefined. The study of second case is more important in the sense 
that the study focuses on deviations between the predictions using a single set of 
coefficients and the act usd bulk flow tnodtl motion as the rotor moves through the 
clearance space in an arbitrary fashion. 

For the purpose of this study, the model of seal represented by a single set of 
coefficients will be identified as lintur tnodtl (e = 0). This model, while valid for a 
small motion about the centered position, may not be accurate for large officenter 
operation of the seal. This off-center motion includes both types of motion described 
earlier. One of the objectives of this study is to identify the magnitude of these 
deviations and examine the effect of these deviations on the overall stability of the 
rotor system and establish limits of effectiveness of using such a model. This task 
is accomplished by solving the bulk flow model seal governing equations directly for 
transient seal forces for any given type of motion, including motion involving large 
eccentricities. Results from this study confirm considerable differences, for large off- 
center operation, between the approximate itnear model (e = 0) and the actual bulk 
flow model 

This approach of solving the governing equations directly for transient seal forces, 
while being the most accurate, may not be practical to be included in a rotor dynamic 
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Figure 7.5 Eccentric Motion: Type 1 



Urn* (iM) 


Figure 7.6 Eccentric Motion: Type 2 
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simulation code, mainly due to the large computing resources required to solve these 
equations at each time step. As a matter of fact, this is the primary reason for using 
approximate models such as the one shown in Eq. (7.1) in rotordynamic simulation 
codes. As an alternative, a general method is developed to model non-linearities 
in an a nn ul ar seal based on dynamic coefficients computed at various static rotor 
operating positions in the seal clearsmce space. This method takes into consideration 
the time history of transient motion, i.e., displacement, velocity and acceleration 
profiles to compute the transient seal forces at any given instant of time. This method 
is extended for approximate displacement, velocity and acceleration profiles to yield 

a practical method that is accurate and easy to implement in a rotordynamic analysis 
code. 

Results from these two methods compare well with those of the actual bulk flow 
model for large eccentric motion. These methods, thoroughly tested for various types 
of transient motion, provide an efficient and practical means for accurate simulation 
of the dynamic eflects of an annular seal for any type of motion. 

The following tasks are accomplished in this study. 

1. Study the efliect of large eccentric motion of the rotor on the dynamic behavior 

®f ® SSME- ATD-liP OTP annular seal using the bulk flow model seal governing 
equations. 

2. Compare the results of the above study with those of the model currently in 
use at NASA/MSFC i.e., model (c = 0). 

3. Develop a method that accurately simulates the dynamics of an «TiTnil«r seal 
for large eccentric motion of the rotor. 

4. Thoroughly test the method for various types of transient motion using bulk 
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flow model results as benclimark. 

5. Compare the results of various models and note their their relative merits and 
deficiencies. 
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CHAPTER Vm 

VARIOUS SEAL MODELS FOR TRANSIENT ANALYSIS 
In tliis chapter, the different models used to study the transient analysis with an 
a nn u lar seal are discussed. These various models differ in the way they compute the 
seal reaction forces for any specified motion of the rotor. The four different models 
used in this study sje explained below. 

1. Bulk Flow Model: This model uses a solution procedure based on the actual 
set of seal governing equations (Eqs. (2.1-2.3)) to compute for transient fluid 
forces at each time step based on a specified motion of the center of the rotor. 
The results of test cases with this model are used as a benchmark to compare 
the other approximate models. 

2. Linear Model (c = 0): This model is based on the linear force-motion model of 
Eq. (7.1) and uses dynamic coefficients computed at zero eccentricity to compute 
the fluid forces. This is the model currently being used at NASA/MSFC for 
SSME turbopump rotordynamic simulations. 

3. New Method-I: This model is based on a new method developed to compute the 

seal forces in a computationally efficient msinner. This method makes 
use of time history of displacement, velocity and acceleration of the rotor to 
compute the seal forces. 

4. New Method-H: This is a simplified extension of method-1 and it assumes 
approximate displacement, velocity and acceleration profiles to compute seal 
forces. 

These various models are shown in flow chart in Figure 8.1. 



Transient Analysis Models 
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Figure 8.1 Various Models for TVansient Analysii 
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In the following sections, each of the above models is discussed in more detail. 


8.1 I^ansient Analysis with Bulk Flow Model 

The linear force-motion model of seal shown in Eq. (7.1) approximates the behavior of 
the bulk flow model governing equations for a small motion of the rotor at the static 
operating point (x,y). The coefficients in this model are obtained by perturbing the 
bulk flow model governing equations given in Eqs (2.1-2.3) at (*,y) and fitting the 
perturbed fluid forces to the linear model. While the linear model is valid only in the 
neighborhood of (x,y), the governing equations are valid at any point in the clearance 
space. To study the deviations between the linear model at a given operating point 
and the actual bulk flow model, this set of governing equations are solved for the 

transient fluid forces directly. The Eqs. (2.1-2.3) are reproduced here for discussion. 
Continuity: 

9{hu) 1 d{hv) dh 

dz R~d0 ^ li ~ ^ 

Axial Momentum: 


pdz - + “al} 

+ + + (V - wy ( 8 . 2 ) 

Circumferential Momentum: 

^ ^ dv^ 

pRdp ~ '•'a + flaa + 

The film thickness in the global coordinate system is given by (Eq. (2.29)) as, 

h(z,/?,<) = ^/(iZ + c)2 - (xstn^ - ycos/?)2 - {xcos0 + ysinfi) - R (8.4) 
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and its derivatives with respect to axial and circumferential coordinates are given by, 


dz 

d0 

di 


iR + c)£_ 


+ cy — {xsinf3 — ycos^Y 
(■R + c)|| — (zstn/3 — ycoafi){zcos0 + ysin/3) 


(8.5) 

( 8 . 6 ) 


^(i? + c)® — {xsinfi — ycos0) 

—{x3in/3 - ycos0){zsin0 - ycos0) ,, ^ . 

j~ s, , . „ ~ {xco»0 + y»in0) (8.7) 

c)^ — {xatn0 — yco80Y 


formula. 



§2 

at 

are approximated using 

a backward difference 

du 


(«(^2) 

-«(<l)) 

(8.6) 

H 


(t2 

-tl) 

dv 


Mi2) 

-^(<l)) 

(8.9) 

di 

ft? 

(t2 

- <l) 

dp 


(piil) 

-P(tl)) 

(8.10) 

di 

ft? 

(t2 



The displacement (x,y) is the displacement of the center of the rotor and (x,y) 

is the velocity of the center of rotor and (i,y) is the acceleration of the center of the 
rotor. 


8.1.1 Transient Seal Forces with Bulk Flow Equations 

The solution procedure to solve the above set of partial differential equations 
is similar to the procedure discussed in Chapter II. At each time step, Eqs. (8.1- 
8.3) are solved to subject to the the boundary conditions given in Eqs. (XX-XX). 
The variables tt(z,/?, t), v(z,/?, f) and p(z,0,t) are the time varying velocities and 
pressure. At each time step, the pressure distribution p(z,0,t) is integrated along 
the length of the seal to compute the two components of the fluid force and 
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^ fluidly* 

- Ftlydd-m{t) = p{z,0,t) COS0 R dfidz (8.11) 

fL #2ir 

- -r/itAJ-vCO = p{z,0,t) sin0 R d0dz (8.12) 

The computation of these transient fluid forces is carried out in a continuous 
fashion from one time step to the next, with the current values of the variables acting 
as initial values for the next time step. Numerical integration with respect to time i 
is implemented using a fourth order Runge-Kutta integrator with adaptive step size. 
The following time step is used for the transient analysis. 

= it (8-13) 

where is the time period of the system. 

8.2 TVansient Analysis with Linear Model (c = 0) 

This is the model usually used to model the dynamic behavior of an annular seal 
(Eq. (7.1)). The dynamic coefficients used in the model refer to those computed at 
the steady state operating position of the rotor. While these linearized coefficients 
can be computed at various eccentricities, it is not possible to decide which set to use 
when the rotor is moving in the clearance space in an arbitrary fashion such as in a 
transient motion. In practice, the set of coefficients computed at zero eccentricity are 
used m the model to compute the seal forces. For example, at NASA/MSFC, SSME 
turbopump simulations use this model. 

8.3 New Method-I 

In this section, a general method is developed to simulate the dynamic behavior of 
an annular seal using djmamic coeffirients computed at various static eccentricities 
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in the clearance space. The motivation for this endeavor is two fold. 

• Firstly, there exists a need for a more accurate model to simulate the dynamic 
behavior of the seal for motion with e > 0.4 as predicted by the bulk flow model 
compared to linear model (c = 0) . Large computing resources required to solve 
the set of governing equations, Eqs. (8.1-8.3), make bulk flow model approach 
impractical for routine rot or dynamic simulations. 

• Secondly, it requires a lot less efibrt to compute dynamic coefficients for a given 
seal umt and for a given set of operating conditions and the methodology for 
this process is well established. 

Consider a single degree of freedom spring-mass-damper system shown in Fig- 
ure 8.2. It is assumed that the stiffness (K), damping (C) and mass (M) vary only with 
displacement * and are independent of velocity and acceleration. This assumption 
follows from the case of an annular seal where the dynamic coefficients are essentially 
functions of eccentricity alone. Various restoring forces in the components of this 
system are considered below. 

8.3.1 Sti&ess Force 

The incremental restoring stiffiiess force A/* in the spring due to an infinitesimal 
extension 6x from * is given by. 


A/t = -K{x)Sx 


(8.14) 


where K is the spring stiffiiess. 
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Figure 8.2 SDOF Spring-Mass-Damper System 

8.3.2 Damping Force 

The incremental restoring damping force A/e in a damper due to an iniinitesimal 
change in velocity 6x is given by, 

A/e = -C{x)Sx (8.15) 

where C is the damping coefficient. 

8.3.3 Inertia Force 

Similarly, the incremental inertia restoring force A/m due to an infinitesimal change 
in acceleration is given by. 

A/m = -M(x)6x (8.16) 



where Af is the inertia or mass coefficient. 
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Figure 8,3 Seal Model for a 2 DOF Vibratioii Model 

8.3.4 Theory 

Consider the 2-DOF model of the seal in Figure 8.3 as represented by the 12 dynamic 
coefficients of the li near force*motion model of Eq. (7.1). at an eccentric position 
(®,y). It is assumed that the stiffiiess, damping and inertia coefficients of this seal 
model arc known at any static eccentricity (x,y) in the seal clearance space. 

Let K„{x^ y), y), k|^(x, y ) and Kyy{x^ y) be the sti&ess coefficients, y ), 

^(**y)i ^w(®jy) the damping coefficients and Afc«(x,y), T7iiry(2fy)i 
T7iy,(x,y), Afyy(x,y) are the inertia coeffidents at eccentridty (x,y). 

An implidt assumption is made regarding the dependence of dynamic coeffidents 
essentially on displacement and they are assumed to be almost independent of ve- 
lodty and acceleration. This assumption will be verified in the following sections in 
comparison with original bulk flow governing equations. 
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The eccentricities or displacements of the center of rotor (x,y) are functions 
of time t and may be specified as (x(t),y(t)). Let and Fy{ti) be the X and 

Y components of the fluid reaction force at any given time f,-. These two force 
components may be considered as a summation of 12 individual component forces 
due to the 12 d3mamic coefficients of Eq. (7.1), as given below. 

-Fy{ti) = /fcv-(i.) + /fcw(*.) + /ev-(i.) + /cyv(«i) + /mv.(<.) + /mw(ti) (8.17) 

Let at tj+j = t,- -I- At, the incremental fluid force components in X and Y directions 
be and AFy{U) respectively. 

/’.(t, + At) = F.(t,) + AF.(tO 

Fy(t, + At) = F„(ti) + AF^{<<) (8.18) 

and the individual components of AF*(tj) and A/^(tj) are given below based on the 
linear-force motion model of Eq. (7.1) 

-AF.(t,) = Af^.{U) + Af,^{ti) + Afa„{ti) + Af^{ti) + Af^{U) + Afr^{U) 

-AFy{U) = A/*^(tO + A/*„(t,)-f A/^(tO + A/.„(t,)-HA/,,,^(t,) + A/„„(t,) 

(8.19) 

The in fini tesimal change in displacement (£zif Sjfi) is pven by, 

= x(ti + At) - x(ti) 


6xx 

hi 


y(ti + At)-y(ti) 


( 8 . 20 ) 
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The mhnitesimal change in velocity is given by, 

Sxi — x{ti + At) — x(ti) 

Syi = y(ti + At) - y(ti) (8.21) 

Similarly, the i nfini tesimal change in acceleration (fxj, 6yi) is given by, 

Sxi — x{ti + At) — x(t,) 

Siii = y{ti + At)-y{ti) (8.22) 

From Eq. (7.1), the incremental fluid forces in terms of their individual components 
are. 


A/s..(tO 



(8.23) 

A/s^(t,) 

= 


(8.24) 

A/fc^(t.) 



(8.25) 

A/iiw(t») 

= 

iifyy(x(tj),y(ij))5yj 

(8.26) 

A/«(t.) 

= 

C„{^{U),y{U))Sxi 

(8.27) 

A/c^(tO 

= 

^ry(®(^0i y(^0)^y» 

(8.28) 


= 


(8.29) 

A/cw(«i) 

= 

Cyy(x(t^),y(tt)) 

(8.30) 

Af^iU) 

= 

Mmm{z{ti),y[ti))6xi 

(8.31) 

A/m»y{t,-) 

= 


(8.32) 

A^p.(t.) 



(8.33) 

A/myv(tt) 


j y(^i))^y* 

(8.34) 
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At time tj+i, after n intervals of At, the incrementiJ fluid force components at each 
time step may be added in the following manner to obtain the total force components. 



~~ 

f=i 

(8.35) 



n 

t=i 

(8.36) 

Ayx(^»+l) 

= 

n 

t=i 

(8.37) 

fkwiU^l) 

= 

n 

(®» » y» ) ^Vi 

t=l 

(8.38) 


= 

n 

t=i 

(8.39) 


= 

n 

Sc^(®My«)^y» 

t=i 

(8.40) 

/ey»(^i+l) 

= 

n 

t=l 

(8.41) 

/ew(^t+l) 

= 

n 

yi)^y» 

t=i 

(8.42) 


= 

n 

t=i 

(8.43) 


= 

n 

, yi )^y* 

«=1 

(8.44) 


= 

n 

isl 

(8.45) 


= 

n 

y»)^y» 

tsl 

(8.46) 


For iniimtesimal quantities At, (fx,^y), (fx, jy) and (5z, fy), the summation may be 
replaced by integration giving the following expressions. 

r"(*) 

= I K„{z,y)dx 
•'0 


( 8 . 47 ) 
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fkxy{t) 


fv(t) 

Jo 

(8.48) 

fkvx{t) 

= 

/•«(») 

kym{x,y)dx 

(8.49) 

/*w(0 

= 

rv(t) 

I ■^w(®»y)^y 

Jo 

(8.50) 

/-.(O 

= 

r*(*) 

J c„{z,y)dx 
Jo 

(8.51) 

/-v(i) 

= 

rv(‘) 

/ <^{x,y)dy 

JO 

(8.52) 

fcy.{t) 

= 

r*(0 

Cyrix, y)dx 

(8.53) 

/cw(0 

= 

rm 

/ Cyy{x,y)dy 

Jo 

(8.54) 

/m.x(0 

= 

r*(‘) 

/ Mrr{x,y)dx 

Jo 

(8.55) 


= 

rv(t) , , 

m^ix,y)dy 

(8.56) 

/cv-(0 

= 

r*(‘) 

(8.57) 

/ew(0 

= 

rii*) 

Myy{x,y)dy 

(8.58) 


These integrals for the case of a single DOF spring*niass*damper system arc 
shown in Figures S.4-8.6. 

8*3.5 Evaluation of Integrals 

Computation of each of the Integrals in Eqs. (8.51-8.61) requires the time history 
of displacement, velocity and acceleration of the rotor center as a function of time 
f. Since each of these curves may have any number of cycles, the following valid 
assumptions are made to reduce summation errors assuming no hystersis loss. 

When the displacement ®(£) is zero the following terms are set to zero and the 
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Figure 8.6 Inertia Force Integral 
integration or summation starts from that instant. 

= 0 ) = 0 

f^{x{t) = 0) = 0 

Similarly, when y{t) is zero the following terms arc set to zero. 

fk^{y{t) = o) = 0 

/fcw(y(0 = 0) = 0 

S im i lar assumptions are made with respect to damping and inertia forces, 
equal to zero the following terms are set to zero. 

/««{x(t) = 0) = 0 

feym{x{t) = 0) = 0 

and at y(t) equal to zero the following terms are set to zero. 


( 8 . 59 ) 

( 8 . 60 ) 
At i(t) 

( 8 . 61 ) 


/«v(y{<) = 0) 


0 
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/ew(y(0 = 0) = 0 (8.62) 

At x(<) equal to zero the following terms are set to zero. 

/m«*(x(t) = 0) = 0 

/myx(®(t) =0) = 0 (8.63) 

and at y(t) equal to zero the foUowing terms are set to zero. 

/m*v(y(t) = 0) = 0 . 

/mw(y(t) = 0) = 0 (8.64) 

The effect of these initializations is to reduce the accumulation of errors as the 
integration is performed along the displacement, velocity and acceleration curves. 
In practice, zero values are never realized and a change in the sign of a variable is 
considered for the above initializations. 


8.3.6 Summation vs. Integration 


In the simulations with this model, the summation approach is used. The time 
step At is made very small so that the summation approaches the integration process. 



(8.65) 


where t„ is the time period of the displacement curve. 


8.3.7 Limitations of Method-I 

The method developed in this section requires the evaluation of 12 integrsJs of 
Eqs. (8.51-8.61). For a general motion, integrating each of these integrals involves two 
highly fluctuating functions. The integrand, which is a dynamic coefficient varying 
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with displacement, and the limit of integration which is either the displacement, 
velocity or acceleration curve. Since closed form solutions are unlikely, one has to 
resort to numerical integration to evaluate these integrals. The accumulation of errors 
over a period of time due to approximate summation or integration schemes as well 
as due to inherent approximate nature of the model may lead to inaccurate results if 
care is no taken to limit these errors. The parameter At is critical for the accuracy 
of this model if summation scheme is used. 

A number of cases using the above method are included in the simulations. 

8.4 New Mcthod-II 

In this section, the method-1 is simplified to make the computation of the Integrals 
in Eqs. (8.51*3.61) easier and more accurate. 

The method- 1 developed in the previous section has the following drawbacks. 

• Difficult to integrate complicated displacement, velocity and acceleration pro- 
files. 

• Accu m ulation of modeling and integration errors as time progresses. 

• Very small time step needed to maintain reasonable accuracy for summation. 

In order to simplify the computation of these Integrals a few assumptions are 
made regarding the displacement, velocity and acceleration profiles. The rrutm as- 
sumption, to be verified, is that the previous time history of the motion has little 
or no effect on the current state of motion for the bulk flow model used. Assuming 
that the above statement is vahd, the actual displacement, velocity and acceleration 
curves may be replaced by approximate curves that are easier to integrate. 
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The following assumptions are made regarding the displacement, velocity and 
acceleration at any given instant of time. 

1. The time history of displacement is neglected and only the current displacement 
is used in the computations. The displacement is assumed to increase &om zero 
to the current value in a linear fashion grsidually, and independent of time. 

2. The time history of velocity prior to the current time step is neglected cmd the 
velocity is assumed to be linear with respect to displacement, i.e., velocity is 
initially zero when the displacement is zero and increases to the current value 
as a li n ea r function of displacement while retaining its direction. 

3. A similar assumption is made about acceleration, i.e., acceleration is zero ini- 
tially and attains its current value in magnitude and direction as a linear func- 
tion of displacement. 

as a function of the time-displacement curve. 

8.4.1 Theory 

As in the previous case, the incremental fluid forces at any given operating point 
are are given by Eq. (7.1). At any given time t*, let the eccentricity or operating 
position of the rotor be given by (e,,^,) or (*i,y,), velocity by or (z‘i,yi) and 

acceleration by (oi,f7j) or (fi,y.). 

The displacement is assumed to incresise from (0,0) to (z,',y,-) linearly. At any point 
along this path, the displacements are given by, 

X = e coa^i 

y = c ain4i (8.66) 
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where, 


e = yjx^ +y* 

e.- = \/*r+yf 

tan^,- = — 

*1 

and velocity as a fonction of displacement is given by, 


(8.67) 


where. 


and acceleration is given as 


X 

= V COS^i 

y 

= i; sin^i 

Vi 

= •/*.■’ + Vi’ 

tanji 

II 

V 

= — e 

function of displacement 


( 8 . 68 ) 


(8.69) 


where. 


X = a cosrfi 

y = a sinrji 


Oi = 

tanrji = 


Oi 
— e 

Ct 


(8.70) 


a 


(8.71) 
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The velocity profile may be rewritten as, 

X = ( — ) X 

y = (-) y (8.72) 

y» 

i.e., the velocity is assumed to increase from (0,0) to the current value (x,-,yi) as a 
linear function of the displacement curve. Similarly, acceleration increases from (0, 0) 
to its current value (xj,yj) as a linear function of displacement. 

X = ( — ) X 

Xi 

y = (-) y (8.73) 

y* 

Using the above assumptions, the Eqs. (8.47—8.58) may be rewritten as, 


fkmmiU) = f K^{x,y) dx (8.74) 

fkmy{ti) = f Kv{x,y)dy (8.75) 

fkvmiU) = ky^{x,y)dx (8.76) 

^w(®jy)<^y (8.77) 

= r C„{x,y){^)dx (8.78) 

fc^iU) = r c^{x,y){^)dy (8.79) 

^0 yi 

fcy^iu) = r Cy,{x,y){^)dx (8.80) 

w 0 

/eyv(^t) = /*" C'w(*»y)(“)^iy (8.81) 

Jo yi 

fm-iu) = J*' M„{x,y){^)dx (8.82) 

fmmviU) — ^>^»»(*»y)(“)<fy 


(8.83) 
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fmymiti) = / my^{x,y){—)dx (8.84) 

JO Xf 

fmyyiU) = ^ Afyy(x, y )( — )dty (8.85) 

•'0 yi 

The integration limits of the above integrals are much simpler as the velocity and 
acceleration profiles arc replaced by equivalent displacement profile which is a simple 
function. These Integrals can be easily computed using any of the various numerical 
integration schemes available. Two methods, one based on Simpson’s rule and other 
based on adaptive quadrature integration are used in this study. 

As compared to the original method, this simplified method has the following 
advantages. 

• Total fluid reaction force is computed at each time step instead of incremental 
summation and this eliminates the problem of accumulation of errors. 

• Integrals are easier to compute. 

• The method is valid for any type of motion. 
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CHAPTER IX 
IMPLEMENTATION 

SimulationE involving various types of transient motion are used to study the 
various approaches discussed in the previous chapter. The model used is a modified 
Jeffcott rotor and is shown in Figure 9.1. It consists of a rotor floating in an *TiTin1»r 
seal and is released &om rest at time f = 0 . Various types of known varying loads 
are applied to the rotor and the resulting motion is studied using the following four 
different models. 

1. Bulk Flow Model: This simulation is done by solving the actual set of seal 
governing equations for transient fluid forces at each time step. The results 
from this study are used as a benchmark to compare the results of the other 
models. 

2. L inear Model (c = 0): This simulation is done using the d}mamic coefficients 
computed at zero eccentricity, to compute the fluid forces. This is the model 
currently in use at NASA/MSFC. 

3. New Method-I: This simulation is done using the method described in section 

8 . 1 . 

4. New Method-II: This is the simplified extension of method-1 and it assumes 
approximate displacement, velocity and acceleration profiles to compute seal 
forces. 
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9.1 Simulation Model 

The equations of motion for the system used for simulation are given by, 

= Fm(t) + 

”»y(0 = ^v(i) + Ffl^.y(t) (9.1) 

where *(<) and y{t) are the rotor displacements from centered position and m is the 
mass of the rotor. The terms Fe(t) and Fy{t) are the two components of the applied 
external load, and Ffiy^^,[t) and Fji^^y{t) are the fluid reaction forces computed 
using one of the four methods mentioned above. 

This set of second order differential equations are reduced to a set of first order 
ordinary differential equations which are then integrated using a fourth order Runge- 
Kutta integrator with adaptive step size. 


9i(0 

= *(<) 


92(<) 

= 


9s(0 

= y{t) 


94(0 

= y(0 

(9.2) 


9i(0 = 92(0 

mgtit) = F,{t) + 

^(0 = 94(0 

mq^{t) = Fy(t) + Ffi^.y{t) (9.3) 

The initial conditions are given below. 


qi(t = 0) = 0 



0 
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q2{t = 0 ) 
qz{t = 0) 
g<{t = 0) 


0 

0 


(9.4) 


At each time step in the above integration process, a routine which computes 
the fluid force is called. The input to the routine is the current displacement, velocity 
and acceleration of the rotor center. The output is the the X and Y components of 
the fluid force. The mass of the rotor used in simulations is 45.4 kg (1001b). 


9.2 Time Step for TVansient Analysis 


In the simulations carried out with various models, the following time step is 
used. 

At=|| (9.5) 

For the method - the following time step is used. 


At = 


_b_ 

500 


About 15-20 cycles of motion is studied for simulation. 


(9.6) 


9.3 Fluid Inertia Coefficients 

The model used to estimate the seal forces has four inertia coefficients 
TOyg and Myy to account for the fluid inertia forces. Of these m^y and are almost 
nero. Strictly speaking, the coefficients and Myy have to be added to the rotor 
mass in the above simulation model. These fluid inertia forces are relatively small 
compared to the other terms even at high frequencies. To simplify the computations 
these inertia coefficients are retained in the linear model and the acceleration at the 
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previouB time step is used to compute these forces. This approach is checked by 
indudiug the fluid inertia coefficients in the rotor mass and comparing the results. 
The results indicate practically no change in the results even at very hi gh frequencies. 

9.4 Computation of Fluid Forces 

In the s imu lations with bulk flow model , the fluid forces are computed by solving 
the set of equations Eqs. (8. 1-8.3) for the pressure distribution p(z,)0,<) which is 
then integrated along the length of the seal to obtain the transient fluid forces. This 
process is implemented in a continuous fashion from one time step to the next. 

For linear model (e = 0) , the set of coefficients for zero eccentricity are used for 
computing the seal forces. These coefficients are obtained from seal code TAMUSEAL- 
III. 

Fot new method-1 and new method-2 , the following procedure is implemented. 
These two methods assume the availability of the 12 d}mamic coefficients as contin- 
uous functions of the displacement (z,y). Let these functions be specified by gimi 
9k*yy 9hyt ctc., as shown below. 

5fc*v(®>y) = ^5iev(®jy) 

yfcv«(®>y) = ^(*iy) 

y*w(®*y) ~ 

9cm,{x,y) = C„(x,y) 

9cmv{x,y) = c^{z,y) 

9eym{x,y) = Cyg{x,y) 
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5cvv(®>y) 

— C|q,(x,y) 

y) 

= M„(z,y) 

y) 

= 

9myx{^^y) 

= nv(*,y) 

ymw(®^ y) 

= -Ai^(z,y) 


(9.8) 


(9.9) 


These carves plotted as a function of eccentricity e are smooth carves with grada- 
ally varying slopes. They have no discontinoities or abrapt changes in function valnes. 
This property of continaoas and smooth variation as a function of eccentricity enables 
these carves to be easily fitted with cubic splines using only a few sets of data. These 
splines can then be interpolated to compute coefficients at any given eccentricity. 

A number of sets of dynamic coefficients at various eccentricity ratios starting 
from 0 to 0.8 are computed using the seal code TAMUSEAL-III . About 10-12 sets of 
dynamic coefficients at 0.05-0.1 eccentricity ratio increments are sufficient for accurate 
interpolation of these coefficients for intermediate values. The increment is made 
smaller for the higher eccentricity region as the coefficients vary more rapidly in that 
region. The data in Table 9.1 is the data used for all the simulations cases. 

These sets of coefficients are obtidned at an eccentridty angle of zero, i.e., along 
the X-axis. Using a transformation, it is possible to compute these coefficients at any 
eccentricity (n,y) using the eccentricity angle <f>. This transformation method has 
been discussed in detml in section 5.4. The same procedure is employed to compute 
d]mamic coefficients in the global coordinate system as required in the simulation 
model of £q. 9.1. 
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9.5 Splines of Coefficients 

About 12 sets of dynamic coefficients computed at various eccentricities and at zero 
eccentricity angle are fitted with cubic splines. This enables the computation of 
the d3ruamic coefficients in a pseudo-continuous fashion. For the case of a distorted 
seal or a seal with a non-circular cross section a 2-D table will be required. The seal 
parameters for seal unit 3-02 as given in Appendix E. A summary of these coefficients 
at various eccentricities is given in Table 9.1. 

9.6 Transient Analysis Simulation Code: TRANSEAL 

The four different models described earlier are implemented in the transient analysis 
simulation code TRANSEAL. The input to this code is the d 3 mamic coefficients at 
various eccentricities, seal parameters and time step At. The output consists of 
the fluid forces, the 12 individual force components, displacements, velocities and 
accelerations at each time step. 



Table 9.1 Table of Dynamic Coefficients for Seal Unit 3-02 
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CHAPTER X 
RESULTS 

To study each of the models described in the earlier chapters, various types of known 
varying loads are applied to a rotor in an unmilAr seal and the resulting transient 
motion is studied. The loads applied to the rotor are divided into the following six 
categories. 

1. Gradually applied loads (ramp function) 

2. Harmonic loads (sinusoidal function) 

3. High frequency loads (sinusoidal function) 

4. Suddenly applied loads (step function) 

5. Impulse or shock loads (impulse function) 

6. Combination of the above loads. 

10.1 Gradually Applied Loads (Ramp Function) 

In this test case, a scries of loads 1780N (4001b), 5340N (12001b) and 8900N (20001b) 
arc applied to the rotor in a vertically downward direction in a gradual manner. 
The loads follow a ramp function while increasing from a zero load at t = 0 to the 
maximum load at t = 0.05s. In practice, this tjrpe of load results from side loads. 
Each of these cases are repeated for the four models and the results are plotted as 
comparison plots. For esich case two plots are shown: the first plot includes the 
time-displacement curve y(f) while the second plot shows the y-component of the 
computed seal force Fy(f) as a function of time. 
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The steady-state values of displacement y marked on the plots are obtained 
separately from a seal code based on the dynamic analysis developed earlier. This 
seal code solves the set of steady-state bulk flow seal governing equations to compute 
these values. These expected steady state displacement values are included for the 
purpose of comparison. 

The comparison plots for this case are shown in Figures 10.1-10.3. 

The external forcing function for this case is given below. 

= -(20t) Fe (0 < f < 0.05s) 

= -Fc (< > 0.05s) 

= 0 ( 10 . 1 ) 

where Fe is the constant load. 

The results for all the three approximate methods show good agreement with 
the bulk flow model for smaller load 1780N (4001b). For higher loads 5340N (12001b) 
and 8900N (20001b), as the eccentricity ratio exceeds 0.4, the 1iti»»AT model (e = 0) 
starts deviating from the actual model and this deviation increases as the eccentricity 
increases. The displacement plot in Figure 10.3 for the 8900N (20001b) case clearly 
shows the difierence between these two models. Also, note the almost exact matching 
of the results of the two new methods with those of the bulk flow model for all cases. 

The plots for the fluid forces in Figures. 10.1-10.3 are all identical because the 
rise time U is much larger than the time period of the system. Therefore, the seal 
reaction forces track the external load exactly in magnitude and without any phase 
shift. 

The following important conclusions can be drawn from these results. 

1. The transient motion using bulk flow model converges to the expected steady 
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Figure 10.1 Gradually AppUed Load, 1780 N (400 lb), Disp. (y), Seal Force (Fy) 
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Figure 10.2 Gradually Applied Load, 5340 N (1200 lb), Disp. (y), Seal Force (Fy) 
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Figure 10.3 Gradually Applied Load, 8900 N (2000 lb), Disp. (y), Seal Force (Fy) 
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state values exactly, thus verifying the transient analysis solution procedure. 
The steady state values used for comp^lrison are obtained independently using 
a separate solution process. 

2. Linear model (c = 0) 3 delds good results for small displacements. At higher 
eccentricities (beyond 0.3— 0.4 eccentricity ratio) there is a marked difference 
between this model and the actual model (Figure 10.3). This example illustrates 
the need for a better model than linear model (c = 0) for motion with large 
eccentricities. 

3. The results of new method-1 and new method-2 exactly match with those of 
the bulk flow model. 

4. On a more important note, unrelated to the transient analysis, these results 
vahdate the solution procedure of the dynamic analysis developed in earlier 
chapters, zeroth order as well as first order solutions. The results of this tran- 
sient analysis establishes the equivalence between the bulk flow model based 
motion and the corresponding motion described by the linear model of Eq.7.1 
based on dynamic coefficients extracted from the first order solution. Generally, 
it is taken for granted that these two solution match for small motion. These 
results demonstrate, for the first time for seals, this equivalence. This exercise 
also may be used as a check case for the first order solution to verify that the 
dynamic coefficients extracted from the first order solution are indeed the cor- 
rect coefficients. In other words, this study may be used as a check case for for 
the solution procedure, zeroth and first order solutions. 
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10.1.1 Steady State Seal Forces vs. Spring Forces 

In a seal analysis, steady state seal forces are computed eks function of eccentricity by 
integrating the zeroth order pressure distribution along the length of the seal. The 
two componenets of the seal force are used to estimate the load bearing capacity of 
the seal. In the context of the new methods 1 and 2, these forces may be thought 
of as forces due to spring stiShesses as damping and inertia forces are non-existent 
when the rotor is in a steady state operating position. The component forces /■ , , 
/k*v» /*vx and f}^ in the new method-1 and new method-2 represent the forces due 
to direct and cross coupled stiffnesses. 

Seal forces F, and Fy shown in Figure 10.4 Seal forces and F„ai-y shown 

in Figure 8.9 are computed using the seal code as a function of eccentricity along the 
X axis. If the new analysis developed is accurate, these forces should be the s£une as 
the spring forces given below. 

Seal forces F, and Fy shown in Figure 10.4 are computed using the seal code 
TAMUSEAL-III as a function of eccentricity along the X axis. If the new analysis 

developed is accurate, these forces should be the same as the spring forces given 
below. 

= /jtcc + fkmy (10.2) 

fkym "t" fkyy (10.3) 

The plot in Figure 10.4 shows the steady state forces computed using both these 
approaches. The results match exactly further validating the two new methods. 
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Figure 10.4 Steady State Seal Forces for Seal Unit 3-02, from Seal Code 
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Figure 10.5 Steady State Seal Forces vs. Spring Forces, New Methods- 1,2 
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10.2 Harmonic Loads (Simisoidal Function) 

In this test case, a series of loads at a frequency of lOOHz are applied to the 
rotor. THs type of forcing function is typical of the loading on the rotor due to an 
unbalance. The external forcing function used for this test case is given below. 

■Fy(t) = F, sin{2itft) 

^-(0 = 0 (10.4) 

where F, is the amplitude and f is the circular frequency in Hz. 

The loads used are 1780N (400)lb and 5340N (12001b) at a frequency of 100 Hz. 
The plots shown in Figures 10.6-10.9 refer to this case. There is excellent agreement 
between the bulk flow model and the all the three approximate models for 1780N 
(4001b) and 3560N (8001b) cases. For the 5340N (12001b) case shown in Figure 10.8, 
the linear model (e = 0) starts to deviate from the actual model. However, the new 
methods gives exact results. 

An important conclusion that can be drawn from this case is that the the linear 
model (c = 0) is accurate even for relatively large displacements provided the rotor 
vibrates about the same static operating point and the set of dynamic coefficients used 
correspond to that operating point. This is interesting since the basic assumption of 

the linear model (e = 0) is that it is valid only for a small motion about the operating 
point. 

For the next case shown in Figure 10.9, the oscillating rotor is forced to move to 
an eccentric position by the application of the following forcing function. 

Fy{t) = (20t) Fc + F, sin{2irft) 0 <t < 0.05s 

= ■P’e + F, sin{2irft) t > 0.05s 
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^*(^) = 0 (10.5) 

where Fg is a constant load. 

This example (Figure 10.9) clearly shows the effect of the movement of the rotor 
operating position away from the centered position as predicted by the bulk flow 
model vis-a-vis the same motion predicted by linear model (c = 0) . There is an 

appreciable difference between these two models. The two new methods give almost 
exact results. 

The following conclusions may be drawn based on these test cases. 

example case in Figure 10.8 refers to the large eccentric motion type-1 
mentioned in Chapter VII. The case in Figure 10.9 refers to large eccentric 
motion type-2. 

2. An annular seal operates almost like a linear element for steady state motion 
about a given static point. 

3. The linear model (e = 0) gives good results even for relatively large rotor 
displacements (e > 0.4). For example, if the rotor is whirling about a static 
operating point, this model yields acceptable results even for large displacements 
if the set of dynamic coefficients corresponding to this operating point are used. 

4. However, the above statement is not true in the cstse of the transient motion 
shown in shown in Figure 10.9 where there is actual movement of the center of 
rotor to am eccentric position. 

5. The two new methods agree exactly with the actual model both for steady state 
as well as trsinsient operation. 
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Figure 10.6 Harmonic Load, 1780N (4001b) at lOOHz, Disp. (y), Seal Force (Fy) 
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= I „ad 3560N (8001b) at lOOHz, Disp. (y), Force (Fy) 

Figure 10.7 Harmonic Load, vouvy , 
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Figure 10.8 Harmonic Load, 5340N (12001b) at lOOHz, Disp. (y), stal Force (Fy) 
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10.3 High Frequency Loads (Sinusoidal Function) 

The rotor of the SSME turbopump runs at about 25000 rpm, i.e., about 417 Hz. Any 
unbalance forces in the rotor system wiU vary as a function of this frequency. To study 
the effect of such unbalance forces at very high frequencies, different harmonic loads 
at 500 Hz and lOOOHz are applied to the rotor and the transient motion is observed. 
The loads applied are 2225N (5001b) and 15570N (35001b) at 500Hz, 17800N (40001b) 
and 3560N (80001b) at lOOOHz. 

The forcing function is similar to the previous case. 


■F’v(t) = F, sin{2vft) 

•^*(0 = 0 ( 10 . 6 ) 

The results are shown in Figures 10.10-10.12 for the 500Hz case and in Fig- 
ures 10.13-10.14 for lOOOHz case. Again, as in the case of harmonic loads at 100 
Hz, there is good correspondence between the bulk flow model and the other three 
methods at both 500 Hz and 1000 Hz cases. Also, these results confirm that the two 
new methods work &t &ny frequency. 

To look at the motion away from the centered position, the rotor is forced to 
an eccentnc position and the results in Figure 10.12 once again clearly show the 
difference the linear model (c = 0) and the actual model. Also, the results of new 
method-1 and new method-2 exactly match those of the actual model. 

The external forcing function for this case is given below. 


^y{t) 

F,{t) 


(20t) Fc -h F, sin{2irft) (0 < t < 0.05s) 
Fe -I- F, sin{2irft) {t > 0.05s) 


0 


(10.7) 
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10.3.1 Fluid Inertia Forces at Higli Frequencies 

The plots shown in Figures 10.10-10.14 may be studied to observe the effect of 
fluid inertia forces at very high frequencies. The results of the three approximate 
methods based on linear model compare very well with the bulk flow model results. 
These results show that the effects of fluid inertia effects at very high frequencies do 
not show any significant difference from the bulk flow model. 

10.4 Suddenly Applied Loads (Step Function) 

In this type of loading, the load is applied suddenly similar to a step function. 
The resulting motion after the transient state settles to the gradually applied load 
case discussed earlier. The loads applied are 1780N (4001b), 3560N (8001b) and 5340N 
(12001b). This motion is almost similar to an impulse or shock load as the rise time 
tr is much smaller than the time period Tn of the system. THs is a good test case 
for the two new methods because of the rapidly varying motion. The results for this 
case are shown in Figures (10.15-10.17). 

The linear model (c = 0) for 5340N (12001b) case in Fig. 10.17 shows that the 
motion not only has a different overshoot but is also quite a bit out of phase and it 
settles to the wrong steady state value. The new method-2 has some overshoot, but 
maintains the phase and settles to the expected steady state value. 

The forcing function for this case is, 

Fy{t) =0 ( i = 0) 

= -F^ (« > 0 ) 

= 0 ( 10 . 8 ) 
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Figure 10.10 High Frequency Load, 2225N (5001b) at 500Hz, Disp. (y) Seal Force 
(Fy) 
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Figure 10.11 High Frequency Load, 15570N (35001b) at 500Hz, Disp. (y). Force (Fy) 
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Figure 10.12 High Frequency Load (eccentric), 4450N at 500Hz, Disp. (y), Force (Fy) 
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Figure 10.13 High Frequency Load, 17800N (40001b) at lOOOHz Disp. (y), Force (Fy) 
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Figure 10.14 High Frequency Load, 35600N (80001b) at lOOOHz, Disp. (y), Force (Fy) 
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Figure 10.15 Suddenly Applied Load, 1780N (4001b), Disp. (y), seal Force (Fy) 
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Figure 10.17 Suddenly Applied Load, 5340N (120( 
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10.5 Impulse Loads (Impulse Function) 

This form of excitation occurs over a very short period of time. Let U be the time 
of force duration and t„ the time period of the system. If t j «: then the applied 

force is classified as an impulse or shock load. Impulse is a force applied over a short 
period of time and is defined as, 

1 = fj‘f{t)dt (10.9) 

where /(<) is the forcing function and I is the impulse. The units of impulse are N-s 
or N-ms. As long as t j the form of f[i) is not important. 

The effect of an impulse on a mass-spring system is to give the mass an initial 
velocity given by, 

m = ^ ( 10 . 10 ) 

and an initial displacement of zero, 

®( 0 ) = 0 ( 10 . 11 ) 

where m is mass of the system. 

In this example case, impulse loads of 2225N-ms (5001b-ms), 4450N-ms (10001b- 
ms) and (18001b-ms) are applied at time t = 0.025s and the motion studied. An im- 
pulse of 2250N-ms (500 lb-ms) refers to an impulse function of a load of 2225N(5001b) 
actmg over a period of 1 ms. Of all the types of loading considered so far, this type of 
load results in maximum overshoot and generally a very rapid varying motion. The 
results for this case are shown in Figures. 10.18-10.21. 

For motion at an eccentricity, the following forcing function is used, with Fe = 
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15001b and / = 5001b-ms. 

F^(0 = (40f)Fc (0<<< 0.025a) 

= Fc + IS{t - 0.50) (t > 0.025a) 

^-(0 = 0 ( 10 . 12 ) 

It is interesting to note that the trsinsient motion predicted by model 

0) in Figure 10.21 is almost 180® out of phase with the actual motion. 

10.6 Combination Loads 

For the final simulation exercise, all the different loads considered in the previous 
simulations are applied simultaneously. The forcing function for this case is given 
below. 

F^(t) = (40t)Fe -I- F,ain(27r/t) (0 < t < 0.025a) 

= Fe + F,atn(27r/i) + IS{t - 0.03) + IS{t - 0.038) (0.025 < t < 0.05a) 
F.(t) = 0 (10.13) 

where 

Fe = 1300 1b 
F, = 1000 lb 
/ = 500 Hz 

I = 500 lb-ms 


and the impulse is applied at t = 0.03s and 0.038s. 

The results for this case are shown in Figure 10.22. 
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10.7 Comparison between Bulk Flow Model and Linear Model (e = 0) 

In Chapter VII, large eccentric motion of the rotor has been classified into the follow- 
ing two types. 

1. Motion with a large amplitude about a static operating point, such as oscillation 
about the centered position. 

2. A similar motion, but with center of the rotor moving in an arbitrary fashion 
away from the center of the seal. 

The results from this study confirm the following general observations for large 
eccentric motion. 

1. Lmear model (e = 0) gives good results as long the the rotor vibrates about a 
static operating point, in this case the centered position, and the set of dynamic 
coeffiaents used in the model correspond to this point. For this type of motion, 
this model is accurate even for relatively large displacements (c > 0.4) and the 
seal behaves like a linear element. 

2. As this operating point moves away from the centered position, the results 
start deviating from the actual model and these deviations get bigger as the 
eccentricity increases, i.e. for e > 0.4. 

10.8 Comparison between Bulk Flow Model and New Method 

The new method-1 and new method-2 give almost similar results for motion with 
smoothly varying displacement, velocity and acceleration profiles. For this type of 
motion the summation technique of new method-1 works well. However, for motion 
with rapidly varying displacement and velocity profiles with abrupt changes in func- 
tion values (sharp peaks and valleys) the summation or integration process encounters 
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problems in some cases. An example of this case is given in Figure 10.23. In this 
simulation, an impulse of 18001b-ms is applied to the rotor and due to the sudden 
change in velocity and acceleration i.e., discontinuities in the motion, the results are 
not accurate. Also, for new method-1 the time step At is critical for accurate evalua- 
tion of the integrals in Eqs. (8.47-8.58) using summation method. The effect of time 
step on the accuracy of this method is shown in Figure 10.24. The first case with a 

fine time step gives exact results. The second case with a coarse time step results in 
an erroneous simulation. 

However, new method-2 is impervious to such problems as integration is carried 

out only with respect to displacement and total fluid force is computed at each time 
step. 

The following observations may be made about these two approaches. 

1. The results show virtually no difference between new method-1 and new method- 
2 , indicating that past motion has practically no effect on the current state of 
motion for the bulk flow model used. 

2. These two methods vastly simplify the computational procedure compared to 
the bulk flow model 

10.9 Conclusions 

This work examines the differences between the linear model (e = 0) and the 
actual bulk flow model for large eccentric motion of the rotor. This study confirms 
considerable deviations between the two modek as the rotor is moved to a large 
eccentric position away from the center of the seal. 

An innovative method to model the seal forces more accurately than the current 
model is developed and tested extensively for various types of transient motion. 
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Figure 10.23 Limitations of Method-I: Discontinuities in Integrals 



Figure 10.24 Limitations of Method-I: Dependence on At 
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The following tasks are accomplished in this study. 

1. Developed transient analysis capability using the original bulk flow governing 
equations. 

2. Established the vahdity of the transient analysis procedure by comparing the 
transient analysis results with the steady state results values obtained form 
TAMUSEAI^III . 

3. Established equivalence between linearized coefficients based transient motion 
and the same motion as predicted by the original governing equations. 

4. The linear model (c = 0) gives good results as long as the rotor vibrates about 
the centered position or a fixed static operating point (with corresponding set of 
dynamic coefficients). This model is vahd even for relatively large displacements 
about this point. For this type of motion, the seal almost acts like a linear 
element. 

5. However, if the rotor operating position moves away from the centered position 
such as in a transient motion, the results show appreciable differences between 
the linear model currently in use and the actual bulk flow model for eccentricity 
ratios above 0.4, the point at which the linear model (6 = 0) starts deviating 

from the actual model, and there exists a need for a more accurate model in 
this region. 

6. Developed a new method and tested it for various types of transient motion. 
This method is vahd for any type of motion including motion at large eccen- 
tricities. 
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7. For the bulk flow model , fluid inertia forces are not significant even at very 
high frequencies. 
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CHAPTER XI 
SEAL CODES 

To implement the analysis developed in earlier chapters, a series of codes were written 
as a part of this research. All the following codes are currently being used at NASA 
Marshall Space Flight Center as primary tools of seal analysis and design of the 

interstage seals of SSMl^ATD-HPOTP. Brief details of these various codes are given 
below. 

11.1 Tamuseal-I 

• Concentric seal analysis. 

• Original Nelson and Nguyen model. 

• Straight, tapered and axially varying profiles. 

• Moody’s and Hirs’ friction models. 

• Constant properties. 

• Very good agreement with check cases. 

• Runs on UNIX based workstation, VAX. 

11.2 Tamuseal-II 

• Eccentric seal analysis 

• Constant properties. 

• Improved dynamic analysis. 



• Moody’s and Hirs’ friction model. 


• Arbitrary profile, axial and circumferential. 

• Distorted seal profile analysis. 

• External load and eccentricity based analysis. 

• Efficient and better mathematical algorithms. 

• Good agreement with check cases. 

• Runs on UNDC based workstation, VAX, and CRAY. 

11.3 Tamuseal-III 

• Eccentric seal analysis. 

• Constant properties. 

• Variable properties from NIST12 (LOX and LH2) 

• Improved dynamic analysis. 

• Moody’s and Hirs’ friction model. 

• Arbitrary profile, axial and circumferential. 

• Distorted profile analysis. 

• External load and eccentricity based analysis. 

• Efficient mathematical algorithms. 

• Good agreement with check cases. 


• Runs on UNIX based workstation, VAX, and CRAY. 
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11.4 Tamuscal-IV 

• Eccentric seal analysis. 

• Constant properties. 

• Variable properties. 

• Tbermal Effects with variable fluid properties. 

• Improved d3Tiamic analysis. 

• Moody’s and Hirs’ friction model. 

• Arbitrary profile, axial and circumferential. 

• Distorted profile analysis. 

• External load and eccentricity based analysis. 

• Eflicient and better mathematical algorithms. 

• Very good agreement with check cases. 

• Runs on UNIX based workstation, VAX, and CRAY. 

11.5 Transeal 

This is the code developed for the study of transient analysis with an annular seal. It 
has the ability to do the following simulations. Given time dependent displacement, 
velocity and acceleration of the center of the rotor, the code computes the transient 
seal forces as a function of time using one of the following four models. 

• Original bulk flow governing equations. 



Linear seal model (e _ 0) with dynamic coefficients computed at zero eccentric- 
ity. 


Using New Method-I. 


Using New Method-II. 
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CHAPTER XII 
CONCLUSIONS 

In this work, a new dynamic analysis for liquid annular seals with arbitrary profile 
is developed based on a method originally proposed by Nelson and Nguyen. The 
foUowing modifications are made to improve the original method. 

1. Improved zeroth order solution based on cubic splines versus FFT method of 
Nelson and Nguyen. The improved method shows better convergence at higher 
eccentricities and yields solution for cases where the original method had failed. 

2. A more exact first order solution based on continuous interpolation of first order 
variables. A new set of equations are derived for dynamic analysis. 

The analysis developed is extended for cryogenic seals for the following models. 

1. Constsmt fiuid properties. 

2. Variable fluid properties, 

3. Thermal efiects (energy equation) with variable fluid properties (concentric 
case). 

A unified solution procedure that is vaKd for both Moody’s fidction model and 
Hirs’ friction model is developed. Dynamic coefficients based on external load speci- 
fication are introduced for seals for the first time. The analysis can be used to model 
seals which support a pre-load. 

Arbitrary profile seals are discussed with reference to an elliptical seal. Unique 
differences that exist between regular straight or tapered circular cross-section seals 
and arbitrary profile seals with a circumferentially varying clearance are analyzed. 
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A study on the effect of orientation of minimuTn film thickness system on computed 
dynamic coefficients is conducted. A general distorted interstage seal of SSME-ATD- 
HPOTP is analyzed as an arbitrary profile seal emd the results are compared to an 
an average clearances profile seal. 

The predictions of current analysis are compared with a number of experimental 
and theoretical cases from seal literature. Based on the comparative studies the 
following conclusions are drawn. 

1 . Good comparison with original Nelson and Nguyen method (1988a, 1988b). Bet- 
ter results for cases where the original method had failed. 

2. Good comparison with Childs’ (1985) Hir’s and Moody’s friction model based 
analyses. 

3. Good comparison with Scharrer and Nelson (1990), except for a discrepancy 
noted in their results. 

4. The deviation between current analysis (variable properties) and San Andres’ 
analysis increase with eccentricity. These differences vary from case to case. 

5. The difference between current analysis (thermal effects with variable proper- 
ties) and San Andres agree well for the concentric case, except for cross coupled 
stifi[he8B. 

The work on transient analysis with an a nnular seal examined the differences 
between linearized coefficients based motion and the actual motion based on bulk 
flow governing equations. This study confirms considerable deviations between the 
two models as the rotor is forced to an eccentric position of c > 0.4. 

Based on this study an equivalence is established between linearized coefficients 
based motion and the same motion based on bulk flow governing equations for small 
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motion. This study may be used to compare the accuracy of solution (dynsunic 
coefficients) of a psirticular dynamic analysis. 

An innovative method to model non-linearities in seal forces is developed. This 
method models seal forces more accurately than the model based on a single set of 
dynamic coefficients. This method is tested extensively for various types of treinsient 
motion. 

The following conclusions are drawn from this study. 

1. Linear model (e = 0) gives good results as long as the rotor vibrates about the 
centered position or a fixed static operating point (with corresponding set of 
d3mamic coefficients). This model is valid even for relatively large displacements 
about this point. For this type of motion, the seal almost acts like a linear 
element. 


2. However, if the rotor operating position moves away from the centered position 
such as in a transient motion, the results show appreciable differences between 
the above linear model and the actual bulk flow model, for eccentricity ratios 
above 0.4. This is roughly the cutoff point where the linear model starts devi- 
ating from the bulk flow model, and there exists a need for a more accurate 
model in this region. 

12.1 Future Work 

One of the areas of improvement for the bulk flow model is the treatment of entrance 
and exit loss coefficients. Typically, in a seal analysis these coefficients are treated as 
constants. It is an accepted fact that these loss coefficients vary with the inlet and 
exit geometries and Reynolds number. A more realistic model for this loss coefficient 
would be a varying coefficient around the inlet and exit circumferences. 
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The other area of improvement for high pressure seals is the inclusion of eleistic 
deformation of the seal housing. Some preliminary studies (Iwatsubo, 1987) have 
already been done in this area. Future work should combine all these various analyses 
into a combined thermoelastic-hydrodynamic (TEHD) analysis. 
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APPENDIX A 

CONSTANT PROPERTIES MODEL 


Zeroth Order (Ste&dy State) Equations: 


where 


Poho 

0 

0 


Bug 

Bz 


F^ 

0 

PohoUo 

0 


Bvq 

Bz 


F. 

PohoUo 

0 

ho 


ex J 


F 

• m 


(A.l) 




— — Potto 


^^0 _ PqVq dho poho dvo 


F„ = 


dz R d0 
_ PohoVQ dvp hodpo 

R W ~ rW 

- Po{frO^\lul + (vo - tw)* + + vg} 

PohpVo dvo 


R df3 

1 

First Order (Perturbed) Equations: Continuity: 


duo , kodvo dho Idho 

rW^ + rW''^ = 


(A.2) 


(A.3) 


- + (vo - wy + f^^yful + vl) (A.4) 


dho 


dho Vo dho 


dt ^ dz R d0 

fdu 1 


Axial Momentum: 


^^*0 ^ l^d^ ^ hoVoduo 

dz 


Po dz ' R dfi 

Circumferential Momentum: 


ji. t, ^**0 , '^tjpo novoouo , , duo 

+ + Ao-^ + AuUi + AfVi = Ahh\ (A.6) 


fc„ti I AflWo dvo ho dpo , dvo _ 

*•“« + -R-gg (A-?) 
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The coefficient expressions etc., in the above first order equations are defined 

below. 


A. 

+ 

s\ ^ 

o 

II 

2 /^•o 

F 

{B»oUtO + Frt 

,Uro) 

(A.8) 

A„ 

_ ho duo 
R dl3 
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r -^90 

Uo{vo— 
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+ (wo-tw)^} 
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Vo dwo 
R dfi 
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+ HMo) 
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ho dvo 
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2^.0 ^ 
’'"UZ* 
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Qs .05 
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Wo dwo 
R d0 





+ + (vo - w)HroUro} (A. 13) 

(A.14) 


with further definitions for Moody’s and Hirs’ friction factors and their dependent 
terms, /»0j/r0>f^»oj-fro etc., given in Appendix D. 

The first-order governing equations are expressed in terms of the ai(z,/3) and 
bi(z,0) functions as; 


ho dni 
p dz 




hoUo dos 
R dl3 

(A.15) 


ho do2 
p dz 


- houas + {A^~ «o|^)a4 + (A, - = 
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~R~d0 


(A.17) 



249 


da^dh ldh ^ ho dae 

fh-X h -^04 + — — 06 = -C^UfCOSfi TT-;^ 

oz dz Rd(3 R 60 


+ -BtiOs + B^Ob + hoUfQo = —C^BhCOS 0 — ~ 
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«o«o-^ + B„04 — houfos + BvOe = 


ho do] hoVo doe 
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(A.19) 

(A.20) 


ho dh\ 
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, dbe , 9h 1 6h vo 6u 1 6v. „ hodbe 

8. + + 58^*“ = + (& + - Ra 0 

*ou.^ + B.I, + BA + *^6. = 

+ Bui4 — houibg + = —c^BhiinP — 


Bp d/3 B 5/? 
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APPENDIX B 

VARIABLE PROPERTIES MODEL 

Zeroth Order (Steady State) Equations: 
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R d/3 


R dpo d/3 


F, = - 


R d^ R d/3 

- Po{/roy -f [vo - «;)* + + 

Po^oVo dvo 

R ~dfi 

- PoifrO ' ° + (vo - + /ri)^\/ugTt^} 


First Order (Perturbed) Equations: Continuity: 
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Axial Momentum: 


hoVo dui dui ho dpi , dui „ ^ 

+ ^o«o~s~ + — + + -B„ui + B„vi + Bppi 


R dp 


dz 


Po dz 
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= Bshi 


(B.6) 
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Circumfereniial Momentum: 


^^0 , ho dpi , , dvi , dvi _ 

R dp ^ PoR dp ^~dt 

= Cnhi (B.7) 


The coefficient expressions Au,A„ etc., in the above first order equations are 
defined below. 
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(B.27) 


with further definitions for Moody’s and Hirs’ friction factors and their dependent 
terms, ftOtfrOiFto^Fro etc., given in Appendix D. 

The first-order governing equations are expressed in terms of the 0,(2, /3) and 
6,(2, y3) functions as; 
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hoVo db^ ho db2 
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APPENDIX C 

THERMAL EFFECTS WITH VARIABLE PROPERTIES MODEL 


Zeroth Order (Steady State) Equations: 


where 
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The first order equations are given in Appendix B. 
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APPENDIX D 

FRICTION FACTORS 


term 

Moody’s Model 

Hirs’ Model 

u 

S^u = + (10* + 10*^)'«) 

fsOB “ 

frO 

Um = ^[1 + (10"^ + 10*^)'^’] 

frOE = WriCcT 

9to 

9 mm- 12H.0 + 

9mH = 

9r0 


9r0E = -rrirnrlRroY^ 

hfo 


hroH = -m,n,[R.o]”‘* 

hrO 

Kom = ^(10“^ + 

Koh = -nv»v[Rro]™' 

F 

/»om/2 

fton/^ 

Fro 

/rOM/2 

frOEl^ 

G,o 

9mmI2 

9meI^ 

GrO 

9tom/^ 

PrOH/2 

H,o 

/»#oif/2 

^»oh/2 

HrO 

^om/2 

KoeI^ 

F.O 

{ftoM — 9 mm)!^ 

iftOE — 9 me)I^ 

Fro 

(/rOM — 9rOM)l2 

ifrOE — 9 t0e)I^ 

U,o 

y/ul + V^ 


UrO 

yjvl + {vo - «;)* 

^1*0 + («0 - V))^ 

Rto 


^y/vl-^vl 

RrO 

+ («„ - «,)■ 

+ (»o - «-)' 
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APPENDIX E 

CLEARANCE FUNCTIONS FOR ELLIPTICAL SEAL 


The clearance functions for the elliptical seal are given in the next table as; 


straight proSle 

linear proSle 

quadratic pro£le 

fi{z) = Oi 

/l(z) = Oi + OjZ 

/i(z) = Cl + ajz + oaz* 

/j(z) = bi 

/2(z) = fci -f bjZ 

/2(z) = 6i + ftjZ + 63Z* 

o 

II 


^ = 02 + 2osz 

o 

II 


^ = 62 + 26 sz 

ai = R + c. 

ui = R + c, 

ai = /Z + C| 


02 = r(Ce-Ci) 

02 = ^(Ce - 4 C„, + 3 C<) 

03 = ^(Oe 2 Cnj "I" Cj) 

= ii + (1 — ^)C| 

6i = iZ + (1 ^)ct 

hi = R + (1 — h)cj 


Is 

1 

1 

-•M 

II 

h2 = + 3 c.) 

hs = X (1 - h)(Ce - 2 c„i + Ci)) 


Gradients of the clearance function for elliptical seal are given by; 


c(z,/3) = 

{fi{z)cosfif + (/ 2 (z)sm^) 

5c 

fif'icos^p + hftsin^^ 

dz 

y/{fiCos/3)^ + (/2Stn/?)’ 

dc 

ifi - ff)cos/3sinl3 




(E.1) 

(E.2) 

(E.3) 
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APPENDIX F 
EQUATION OF STATE 


The modified Benedict- Webb-Ruben Equation of state is given below. 

P = pRT + 

P^{G{1)T + -h (7(3) (?(4)/T -h G{5)/T^) + 

P^{Gi6)T + G{7) + C?(8)/r -h C?(9)/T*) -i- 

P^(C?(10)r + ^(11) -h C7(12)/r) -h /?®(G(13)) -I- 

/(C?(14)/T + C7(15)/r’) -h p''{G{16)/T) -h 

p\G{n)!T + (?(18)T*) -I- p®(C7(19)/r*) -1- 

A>^(C?(20)/r* -h G(21)/T*) + 

p®(C?(22)/r* + G(23)/r^) -t- 

/(G(24)/r* + G(25)/T*) -i- 

/(G(26)/r* -h G{27)fT*) + 

p“(G(28)/T* + G(29)/r*) -t- 


p'®(G(30)/r’ -i- G(31)/T* -I- G(32)/2^) (F.l) 

The expression for viscosity is given by, 

+ fii{T)p + fii(pjT) (F.2) 

» = i; G.(i)r(‘-')/> (F.3) 

t=l 

/ii(T) = F.(l) + F,(2) {F.(3) - Jn(7’/F.(4))}’ (F.4) 

Mp,T) = £'■<'•’■1 - e®m fF.5) 


F(P.T) = 


£.( 1 ) + E.( 2 )ir(p) + E,(3y‘ + 
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E^H)H{p)IT’ + £,(6)/Vr‘' + 
E,(6)/T + E,{7)H(p)/T 
G{T) = £„( 1 ) + K(2)/r 

H{p) = p°'(p - E^(S))/E4S) 


(F.6) 

(F.7) 

(F.8) 


The variation of density with respect to pressure (isothermal case), ^ (1/||) is 
given as, 

dp 

^ = RT + 

dp 

2p{G{l)T + G(2)T'/* + G(3) + G'(4)/T + C7(5)/T*) + 
3p*(G(6)T + G{7) + G(8)/T + G(9)/r*) + 

V(G(10)T + G(ll) + G(12)/T) + 5p^(G(13)) + 

6p®(G(14)/T + G(15)/T*) + 7p«(G(16)/T) + 

8p^(G(17)/r + G(18)T’) + 9p*(G(19)/r*) + 

{3p2 + 27p^}(G(20)/r* + G(21)/T*) + 

{5p* + 27p®}(G(22)/r* + G{2Z)!T*) + 

{7p* + 27 p*}(G(24)/T* + G(25)/T») + 

{9p» + 27 p^0}(G(26)/T’ + G(27)/2^) + 

{1 Vo + 27 p^2}(G(28)/T’ + G(29)/T*) + 

{13p'2 + 27p'4}(G(30)/r* + G(31)/T’ + G(32)/7^) (F.9) 


The variation of viscosity with respect to pressure (isothermal case), ^ 

is ^ven below. 


= ,i.(r) + + E.(4)i + £.(7)i 


dp 
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dH 

dp 


+ + £45)5^}] 


(F.IO) 

(F.ll) 


Ap\ 

A^i 

A^i 


= ^ 
dp 

= ^ 
dp 

_ ^ 
dp 

= X B^x 


(F.12) 

(F.13) 

(F.14) 

(F.15) 
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LIQUID SEAL DATA 
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The following data is for Childs and Lindsey (1993) 


Seal Data for 

Childs/Lindsey (1993) Cast 

seal length, L 

13.13 mm 

rotor radius, R 

76.20 mm 

c. 

nominal values are given in Chapter VI 

Ce 

for measured values see Lindsey (1993) 

fluid 

water at 54.5® C 

density, p 

985.25 kg/m® 

viscosity p 

0.5268x10-® Pa-s 

pressure drop, Ap 

1.38 Mpa, 2.41 Mpa, 3.45 Mpa 

rotor speed, N 

10200 rpm, 17400 rpm, 24600 rpm 

friction factor 

Moody Model 

relative roughness Cr/2c, 

0.001 (rotor) 

relative roughness e,/2c. 

0.001 (stator) 

pre-swirl ratio 

see table in Lindsey (1993) 

inlet loss, 

0.1 

exit pressure recovery, 

1.0 




The following data is for Childs and Kim (1985) case. 


Seal Data for 

Childs & Kim (1985) Case 

seal length, L 

16.66 mm (0.656 in) 

rotor radius, R 

48.39 mm (1.905 in) 

Ci 

0.224 mm (0.0088 in) 

c. 

0.158 mm (0.0062 in) 

C* 

0.191 mm (0.0075 in) 

fluid 

liquid oxygen 

density, p 

1121.26 kg/m^ (70.0 Ibm/ft’) 

viscosity pi 

1.292x10"'* Pa-s (2.70x10"* Ib-s/ft*) 

pressure drop, Ap 

32.26 MPa (4700 psi) 

rotor speed, N 

23700 rpm 

friction factor 

Hirs Model 

roughness 174 , 

-0.0980, 0.01091 (rotor) 

roughness m„ n. 

-0.0433, 0.03120 (stator) 

pre-swirl ratio 

0.3 

inlet loss, 

0.1 

exit pressure recovery, 

1.0 
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The following data is for Sharrer and Nelson (1990) case. 


Stal Data for 

Sharrer & Nelson (1990) Case 

seal length, L 

22.20 nun (0.874 in) 

rotor radius, R 

42.50 mm (1.673 in) 

Ci 

0.381 nun (0.015 in) 

c. 

0.127 nun (0.005 in) 

C* 

0.127 nun (0.005 in) 

fluid 

liquid oxygen 

density, p 

1124 kg/m* (70.0 Ibm/ft*) 

viscosity p, 

1.34x10-“ Pa-s (2.70x10-* Ib-s/ft*) 

pressure drop, Ap 

44.82 MPa (6500 psi) 

rotor speed, N 

30000 rpm 

friction factor 

Hirs Model 

roughness 

-0.2500, 0.0790 (rotor) 

roughness m,,n. 

-0.1360, 0.0697 (stator) 

pre-swirl ratio 

0.8 

inlet loss, 

0.25 

exit pressure recovery, ^ 

, 1.0 




The following data is for Scharrer and Nunez (1989) rase 


Stal Data for 

Sharrer & Nunez (1989) Case 

seal length, L 

45.70 mm (1.8 in) 

rotor radius, R 

45.50 mm (1.79 in) 

nominal Ci 

0.107 mm (0.0042 in) 

nominal Ce 

0.089 mm (0.0035 in) 

c. 

0.089 mm (0.0035 in) 

fluid 

liquid oxygen 

density, p 

1124 kg/m» (70.0 Ibm/ft’) 

viscosity fi 

1.34x10"^ Pa-s (2.70x10"® Ib-s/ft*) 

pressure drop, Ap 

14.93 MPa (2165 psi) 

rotor speed, N 

37360 rpm 

friction factor 

Hirs Model 

relative roughness m„n„ 

-0.25, 0.079 (rotor) 

relative roughness m„n. 

-0.136,0.0697 (stator) 

pre-swirl ratio 

0.6 

inlet loss, 

0.25 

exit pressure recovery, 

1.0 
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The following data is for Jenssen (1970) case. 


Seal Data for Jenssen (1970) Case 

seal length, L 

50.80 nun (2.0 in) 

rotor radius, R 

24.76 mm (0.975 in) 

Ci 

0.272 mm (0.0107 in) 

Ce 

0.272 mm (0.0107 in) 

€• 

0.272 mm (0.0107 in) 

fluid 

water 

density, p 

1000 kg/m» (62.43 Ibm/ft*) 

viscosity p, 

1.30x10-® Pa-s (2.72x10-® Ib-s/ft®) 

pressure drop, Ap 

0.344, 1.034, 1.724 MPa (50,150,250 psi) 

rotor speed, N 

3000, 5000, 7000 rpm 

friction factor 

Moody Model 

relative roughness e,/2c. 

0.0 (rotor) 

relative roughness e,/2c, ^ 

0.0000748 (stator) 

pre-swirl ratio 

0.3 

inlet loss, 

0.4 

exit pressure recovery, 

1.0 




The following data is for Falco et al. (1986) case. 


Seal Data for Falco et al. (1986) Case 

seal length, L 

40 mm (1.57 in) 

rotor radius, R 

80 mm (3.15 in) 

c» 

0.36 mm (0.0142 in) 

c. 

0.36 mm (0.0142 in) 

C. 

0.36 mm (0.0142 in) 

flmd 

water 

density, p 

1000 kg/m’ (62.43 Ibm/ft’) 

viscosity p 

1.0x10-’ Pa-s (2.10x10-* Ib-s/ft’) 

pressure drop, Ap 

1.0 MPa (145 psi) 

rotor speed, N 

4000 rpm 

friction factor 

Moody Model 

relative roughness e, /2c, 

0.0044 (rotor) 

relative roughness e,/2c. 

0.0083 (stator) 

pre-swirl ratio 

0.3 

inlet loss, 

0.3 

exit pressure recovery, 

1.0 




The following data is for Kanki sind Kawakami (1984) case. 


Stal Data for Kanki & Kawakami (1984) Cast 

seal length, L 

200 mm (7.87 in) 

rotor radius, R 

100 mm (3.94 in) 

c. 

0.50 mm (0.0197 in) 

c« 

0.50 mm (0.0197 in) 

C* 

0.50 mm (0.0197 in) 

fluid 

water 

density, p 

1000 kg/m® (62.43 lbm/ft») 

viscosity fi 

1.0x10“’ Pa^B (2.10x10“’ Ib-s/ft’) 

pressure drop, Ap 

0.98 MPa (142 psi) 

rotor speed, N 

2000 rpm 

friction factor 

Moody Model 

relative roughness e,/2c. 

0.0033 (rotor) 

relative roughness e,/2c. 

0.0033 (stator) 

pre-swirl ratio 

0.0 

inlet loss, 

0.1 

exit pressure recovery, 

1.0 




The following data is for AUsdre et al. (1976) case. 


Seal Data fc 

Allaire et al. (1976) Case 

seal length, L 

40.6 mm (1.60 in) 

rotor radius, R 

39.9 mm (1.57 in) 

Ci 

0.14 mm (0.0055 in) 


0.14 mm (0.0055 in) 

C* 

0.14 TTim (0.0055 in) 

fluid 

liquid hydrogen 

density, p 

57.657 kg/m^ (3.6 Ibm/ft’) 

viscosity p 

7.4396x10-* Pa-s (1.5538xl0~^ Ib-s/ft*) 

pressure drop, Ap 

7.26 MPa (1050 psi) 

rotor speed, N 

23700 rpm 

friction factor 

Moody Model 

relative roughness e,. /2c, 

0.00 (rotor) 

relative roughness e,/2c. 

0.000001 (stator) 

pre-swirl ratio 

0.5 

inlet loss, 

0.1 

exit pressure recovery, 

1.0 
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The following data is for San Andres et al. (1992) 


Stal Paramtttrs for San Andres et al. (1992) 

seal length, L 

45.70 mm (0.656 in) 

rotor radius, R 

45.50 mm (1.905 in) 


0.127 mm (0.00587 in) 

Ce 

0.127 mm (0.00581 in) 


0.127 mm (0.00587 in) 

fluid 

liquid oxygen 

density, p 

variable prop. frx>m NIST12 (MIPROPS) 

viscosity p 

variable prop, from NIST12 (MIPROPS) 

inlet pressure, ir 

18.31 MPa (2621 psi) 

exit pressure, p« 

3.378 Mpa (483 psi) 

rotor speed, N 

37360 rpm 

friction factor 

Moody Model 

relative roughness, e, /2c, 

0.0 (rotor) 

relative roughness, e,/2c. 

0.044 (stator) 

pre-swirl ratio 

0.6 

inlet loss, 

0.25 

exit pressure recovery, 

1.0 
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The following data is for elliptical seal case. 


Seal Parameters for Elliptical Seal 

seal length, L 

16.66 mm (0.656 in) 

rotor radius, R 

48.39 mm (1.905 in) 

c. 

0.069 mm (0.00273 in) 

Ce 

0.099 mm (0.00390 in) 

C* 

0.069 mm (0.00273 in) 

fluid 

liquid oxygen 

density, p 

1041.7 kg/m’ (65.03 Ibm/ft’) 

viscosity p 

1.296x10-“ Pa-s (0.188 xlO"’ Ib-s/ft’) 

pressure drop, Ap 

25.39 MPa (3681 psi) 

rotor speed, N 

22700 rpm 

friction factor 

Moody Model 

relative roughness, Cr/2c. 

0.0 (rotor) 

relative roughness, e,/2c. 

0.03 (stator) 

pre-swirl ratio 

0.2 

inlet loss, 

0.33 

exit pressure recovery, 

1.0 




The foUowing data is for the distorted seal case. 


Seal Parameters for Distorted Seal Unit S-01 

seal length, L 

16.66 nun (0.656 in) 

rotor radius, R 

48.39 mm (1.905 in) 

average c, 

0.149 mm (0.00587 in) 

average c. 

0.148 mm (0.00581 in) 

nominal clearance c* 

0.149 mm (0.00587 in) 

fluid 

liquid oxygen 

density, p 

1041.7 kg/m» (65.03 lbm/ft>) 

viscosity /i 

1.296x10"“ Pa-s (0.188 xlO"* Ib-s/ft’) 

pressure drop Ap 

35.25 MPa (5112 psi) 

rotor speed, N 

25000 rpm 

friction factor 

Moody Model 

relative roughness, e, /2c, 

0.0 (rotor) 

relative roughness, e,/2c. 

0.8518 (stator) 

pre-swirl ratio 

0.2 

inlet loss, 

0.3 

exit pressure recovery, 

1.0 
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The following data is for seal unit 3-02 used in transient simulations. 


Seal Parameters for Seal Unit 3-02 

seal length, L 

16.6 mm (0.656 in) 

rotor radius, R 

45.7 mm (1.80 in) 

Ci 

0.174 mm (0.00687 in) 

Ce 

0.148 mm (0.00581 in) 

Bominal clearance c« 

0.149 mm (0.00581 in) 

fluid 

liquid oxygen 

density, p 

1041.7 kg/m® (65.03 lbm/ft») 

viscosity fi 

1.296x10-^ Pa-s (0.188xl0~* Ib-s/ft*) 

pressure drop Ap 

35.25 MPa (5112 psi) 

rotor speed, N 

25000 rpm 

friction factor 

Moody Model 

relative roughness, e, /2c, 

0.0 (rotor) 

relative roughness, e,/2c. 

0.8518 (stator) 

pre-swirl ratio 

0.2 

inlet loss, 

0.3 

exit pressure recovery, 

1.0 





